remarkRemark exampleExample

# Stability of mixed FEMs for non-selfadjoint indefinite second-order linear
elliptic PDEs

C. Carstensen Humboldt-Universität zu Berlin, Berlin, Germany & Distinguished
Visiting Professor, Department of Mathematics, Indian Institute of Technology
Bombay, Powai, Mumbai 400076, India (cc@math.hu-berlin.de). Neela Nataraj
Department of Mathematics, Indian Institute of Technology Bombay, Powai,
Mumbai-400076, India (neela@math.iitb.ac.in). Amiya K. Pani Department of
Mathematics, Indian Institute of Technology Bombay, Powai, Mumbai-400076,
India (akp@math.iitb.ac.in).

###### Abstract

For a well-posed non-selfadjoint indefinite second-order linear elliptic PDE
with general coefficients $\mathbf{A},\mathbf{b},\gamma$ in $L^{\infty}$ and
symmetric and uniformly positive definite coefficient matrix $\mathbf{A}$,
this paper proves that mixed finite element problems are uniquely solvable and
the discrete solutions are uniformly bounded, whenever the underlying shape-
regular triangulation is sufficiently fine. This applies to the Raviart-Thomas
(RT) and Brezzi-Douglas-Marini (BDM) finite element families of any order and
in any space dimension and leads to the best-approximation estimate in
$H(\operatorname{div})\times L^{2}$ as well as in in $L^{2}\times L^{2}$ up to
oscillations. This generalises earlier contributions for piecewise Lipschitz
continuous coefficients to $L^{\infty}$ coefficients. The compactness argument
of Schatz and Wang for the displacement-oriented problem does not apply
immediately to the mixed formulation in $H(\operatorname{div})\times L^{2}$.
But it allows the uniform approximation of some $L^{2}$ contributions and can
be combined with a recent $L^{2}$ best-approximation result from the medius
analysis. This technique circumvents any regularity assumption and the
application of a Fortin interpolation operator.

Keywords mixed finite element method, stability, non-selfadjoint indefinite,
general linear second-order elliptic PDE, RT and BDM finite elements, best-
approximation, medius analysis

AMS subject classification 65N12, 65N15, 65N30

##  1 Introduction

This section introduces the non-selfadjoint indefinite second-order linear
elliptic PDE and its mixed formulations. A brief review of earlier results is
followed by the assertion of the stability and the best-approximation results.

###  1.1 Non-selfadjoint indefinite second-order linear elliptic PDEs

The strong formulations for second-order elliptic problems with coefficients
${\bf A},$ b, $\gamma$ componentwise in $L^{\infty}(\Omega)$ and $f\in
L^{2}(\Omega)$ read $\mathcal{L}_{j}u_{j}=f$ a.e. in a polyhedral bounded
Lipschitz domain $\Omega\subset\mathbb{R}^{n}$ with homogeneous Dirichlet
boundary condition $u_{j}=0$ on $\partial\Omega$ for $j=1,2$ and any dimension
$n\geq 2$. For all $v\in H^{1}_{0}(\Omega)$, the two differential operators
(referred to as conservative resp. divergence form throughout this paper) read

(1) $\displaystyle\mathcal{L}_{1}v:=-\nabla\cdot(\mathbf{A}\nabla
v+v\,{\mathbf{b}})+\gamma\,v\quad\text{and}\quad\mathcal{L}_{2}v:=-\nabla\cdot(\mathbf{A}\nabla
v)+{\mathbf{b}}\cdot\nabla v+\gamma\,v.$

The assumption on ellipticity means that the $n\times n$ coefficient matrix
$\mathbf{A}(x)$ is symmetric and positive definite with eigenvalues in one
universal compact interval of positive reals for a.e. $x\in\Omega$. This makes
$\mathcal{L}_{1},\mathcal{L}_{2}:H^{1}_{0}(\Omega)\to H^{-1}(\Omega)$ Fredholm
operators of index zero and their weak formulations
$a(v,w):=\langle\mathcal{L}_{1}v,w\rangle_{H^{-1}(\Omega)\times
H^{1}_{0}(\Omega)}=\langle\mathcal{L}_{2}w,v\rangle_{H^{-1}(\Omega)\times
H^{1}_{0}(\Omega)}$, for all $v,w\in H^{1}_{0}(\Omega)$, are dual to each
other in the duality bracket
$\langle\bullet,\bullet\rangle_{H^{-1}(\Omega)\times H^{1}_{0}(\Omega)}$ of
$H^{-1}(\Omega)$, the dual of $H^{1}_{0}(\Omega)$.

Throughout this paper, zero eigenvalues are excluded and the kernel (of one of
these operators) $\mathcal{L}_{j}$ is supposed to be trivial, so that
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ are bijections. It is known from the
theory of bilinear forms in reflexive Banach spaces [2, 3] that this implies
well-posedness and the continuous $\inf$-$\sup$ condition (e.g., when
$H^{1}_{0}(\Omega)$ is endowed with the norm $\|\nabla\bullet\|$)

(2) $0<\alpha:=\inf_{v\in H^{1}_{0}(\Omega)\setminus\\{0\\}}\sup_{v\in
H^{1}_{0}(\Omega)\setminus\\{0\\}}\frac{a(v,w)}{\|\nabla v\|\,\|\nabla w\|}.$

The $\inf$-$\sup$ constant is the same for the original and the dual problem;
$a(v,w)$ could be replaced by $a(w,v)$ with the same $\alpha$. The finite
element error analysis is enormously simplified under additional conditions on
the coefficients that lead to an ellipticity of $a(\bullet,\bullet)$ and allow
an application of the Lax-Milgram lemma [2, 3, 4]. The present situation of a
general non-selfadjoint indefinite second-order linear elliptic PDE avoids any
of those assumptions and examines coefficients in $L^{\infty}$, which satisfy
the following.

Assumption (A). There exist two global constants
$0<\underline{\alpha}\leq\overline{\alpha}<\infty$ such that ${\bf A}\in
L^{\infty}(\Omega;\mathbb{R}^{n\times n})$ satisfies
$\underline{\alpha}\leq\lambda_{1}({\bf A}(x))\leq\cdots\leq\lambda_{n}({\bf
A}(x))\leq\overline{\alpha}$ for the eigenvalues $\lambda_{1}({\bf
A}(x))\leq\dots\leq\lambda_{n}({\bf A}(x))$ of the SPD ${\bf A}(x)$ for a.e.
$x\in\Omega$. The functions $\mathbf{b},\mathbf{b}_{1},\mathbf{b}_{2}\in
L^{\infty}(\Omega;\mathbb{R}^{n})$ and $\gamma\in L^{\infty}(\Omega)$ are
componentwise bounded in the bounded polyhedral Lipschitz domain
$\Omega\subset\mathbb{R}^{n}$.

Given the various applications to porous media and ground-water flow with
rough and oscillating coefficients merely bounded in a well-posed PDE, this
contribution gives an affirmative answer to the fundamental question whether
the mixed finite element method be used (and then is stable and provides best-
approximation property at least for fine triangulations).

###  1.2 Earlier contributions

For conforming finite element discretizations and sufficiently small mesh
sizes, [17] establishes the existence and uniqueness of conforming finite
element solutions under assumption (A). The mixed formulation for the
conservation (resp. divergence) equation $\mathcal{L}_{1}u=f$ (resp.
$\mathcal{L}_{2}u=f$) introduces the flux variable $\sigma=-\mathbf{A}\nabla
u-u\,{\mathbf{b}}$ (resp. $\sigma=-\mathbf{A}\nabla u$) and seeks the solution
$x=(\sigma,u)\in H$ to

(3) $b(x,y)=(f,v)_{L^{2}(\Omega)}\quad\text{for all }y=(\tau,v)\in
H:=H(\operatorname{div},\Omega)\times L^{2}(\Omega)$

with ${\mathbf{b}}_{1}:=\mathbf{A}^{-1}{\mathbf{b}}$, ${\mathbf{b}}_{2}:=0$
(resp. ${\mathbf{b}}_{1}:=0$,
${\mathbf{b}}_{2}:={\mathbf{b}}\cdot\mathbf{A}^{-1}$) and
$(\sigma,\tau)_{\mathbf{A}^{-1}}:=(\mathbf{A}^{-1}\sigma,\tau)_{L^{2}(\Omega)}$
in

(4) $\displaystyle b(x,y)$
$\displaystyle=(\sigma,\tau)_{\mathbf{A}^{-1}}-(u,\operatorname{div}\tau)_{L^{2}(\Omega)}+(v,\operatorname{div}\sigma)_{L^{2}(\Omega)}$
$\displaystyle+(u,\mathbf{b}_{1}\cdot\tau)_{L^{2}(\Omega)}-(v,\mathbf{b}_{2}\cdot\sigma)_{L^{2}(\Omega)}+(\gamma\,u,v)_{L^{2}(\Omega)}.$

The equivalence to the boundary value problems associated with the linear
differential operators in (1) and their well-posedness on the continuous level
can be found in [5, Sect. 2]. This implies the continuous $\inf$-$\sup$
conditions [2, 3]

(5) $0<\beta:=\inf_{x\in H\setminus\\{0\\}}\sup_{y\in
H\setminus\\{0\\}}\frac{b(x,y)}{\|x\|_{H}\|y\|_{H}}=\inf_{y\in
H\setminus\\{0\\}}\sup_{x\in
H\setminus\\{0\\}}\frac{b(x,y)}{\|x\|_{H}\|y\|_{H}}.$

The existence and uniqueness of discrete solutions and optimal $L^{2}$ error
estimates were introduced in [9] for sufficiently fine triangulations in two
and three space dimensions under high regularity assumptions, where the pair
$(\sigma_{h},u_{h})$ is approximated in $RT_{k}(\mathcal{T})\times
P_{k}(\mathcal{T})$ with the Raviart-Thomas (RT) for $2D$ (resp. Raviart-
Thomas-Nedelec for $3D$) finite elements. Global $L^{\infty}$ and global
$L^{2}$ and negative norm estimates for the conservation form were discussed
in [10, 14] for smooth coefficients.

Provided the coefficients $\mathbf{A}$ and ${\mathbf{b}}$ are Lipschitz
continuous, $\gamma$ is piecewise Lipschitz continuous and $H^{2}$ regularity
of the adjoint system, an interesting convergence phenomenon for the BDM
finite element family is clarified in the fairly general framework of [8].

Let $M_{k}(\mathcal{T})$ be any RT or BDM finite element space of degree
$k\in\mathbb{N}_{0}$ and define the discrete space
$V(\mathcal{T}):=M_{k}(\mathcal{T})\times P_{k}(\mathcal{T})\subset H$ based
on a shape-regular triangulation $\mathcal{T}$ with mesh-sizes $\leq\delta$,
written $\mathcal{T}\in\mathbb{T}(\delta)$.

In case $\mathbf{A}$ and $\mathbf{b}$ are globally Lipschitz continuous and
$\gamma$ is piecewise Lipschitz continuous, the convergence results in [8]
also establish stability in the sense

(6) $0<\beta_{0}\leq\inf_{\mathcal{T}\in\mathbb{T}(\delta)}\inf_{x_{h}\in
V(\mathcal{T})\setminus\\{0\\}}\sup_{y_{h}\in
V(\mathcal{T})\setminus\\{0\\}}\frac{b(x_{h},y_{h})}{\|x_{h}\|_{H}\|y_{h}\|_{H}}(=:\beta_{h})$

for some positive $\delta$ and $\beta_{0}$. With extra work and refined
arguments along the lines of [8], but with reduced elliptic regularity and
solution $u_{j}\in H^{1}_{0}(\Omega)\cap H^{1+s}(\Omega)$ to
$\mathcal{L}_{j}u_{j}=f$ for some $s>0$. Those arguments are not valid under
Assumption (A).

Modern trends in the mathematics of mixed finite element schemes include local
stable projections with commuting properties [11, 12, 13]; those techniques do
not seem to allow the proof of discrete stability and best-approximation under
assumption (A).

Piecewise Lipschitz continuous coefficients with regularity in $H^{1+s}$ (for
some positive $s$) lead in [5] to stability for the lowest-order RT FEM. The
equivalence to nonconforming Crouzex-Raviart finite elements holds more
generally [1] and the combination with the arguments from [17] and [5] might
lead to stability results under the assumption (A) for more examples. In
comparison, the methodology of this paper provides stability for any degree
$k$ and any dimension $n$ (RT and BDM merely serve as popular model examples).

###  1.3 Contribution of this paper

Under the Assumption (A) and for any RT or BDM finite element space
$V(\mathcal{T}):=M_{k}(\mathcal{T})\times P_{k}(\mathcal{T})\subset
H:=H(\operatorname{div},\Omega)\times L^{2}(\Omega)$ of degree
$k\in\mathbb{N}_{0}$ [2, 3, 4], the discrete stability (6) is established for
small mesh-sizes, where either ${\mathbf{b}}_{1}:=\mathbf{A}^{-1}{\mathbf{b}}$
and ${\mathbf{b}}_{2}:=0$ or ${\mathbf{b}}_{1}:=0$ and
${\mathbf{b}}_{2}:={\mathbf{b}}\cdot\mathbf{A}^{-1}$ in (4).

######  Theorem 1.1 (discrete stability).

For each (positive) $\beta_{0}<\beta$ with $\beta$ from (5), there exists
$\delta>0$ such that (6) holds.

This theorem implies [2, 3] that the mixed finite element problems for the RT
and the BDM finite element families of any degree $k$ and in any space
dimension $n$ are (i) uniquely solvable, (ii) uniformly bounded in $H$, and
(iii) fullfil quasi-optimal error estimates in the norm of $H$, whenever the
underlying shape-regular triangulation is sufficiently fine.

Theorem 1.1 and the tools of this paper lead to $L^{2}$ best-approximation up
to oscillations.

######  Theorem 1.2 ($L^{2}$ best approximation).

Suppose $\delta>0$ satisfies (6) with $b(\bullet,\bullet)$ defined for general
$\mathbf{b}_{1},\mathbf{b}_{2}\in L^{\infty}(\Omega;\mathbb{R}^{n})$ under
Assumption (A). Assume $\mathcal{T}\in\mathbb{T}(\delta)$ and that
$x:=(\sigma,u)\in H$ (resp. $x_{h}\equiv(\sigma_{h},u_{h})\in
V_{h}:=V(\mathcal{T}):=M_{k}(\mathcal{T})\times P_{k}(\mathcal{T})$) satisfy
$b(x-x_{h},y_{h})=0$ for all $y_{h}\in V_{h}$. Then the following results (a)
and (b) hold.
(a) There exists a positive constant $C_{1}$, which exclusively depends on
$\beta_{0}>0$, the $L^{\infty}$ norms of (all the components of)
$\mathbf{A}^{1/2}\mathbf{b}_{1}$, $\mathbf{A}^{1/2}\mathbf{b}_{2}$, and
$\gamma$, as well as on the shape-regularity of $\mathbb{T}$, such that the
piecewise mesh size $h_{\mathcal{T}}$ in $\mathcal{T}$ and the $L^{2}$
projection $\Pi_{k}$ onto $P_{k}(\mathcal{T})$ satisfy

$C_{1}^{-1}\left(\|\sigma-\sigma_{h}\|_{\mathbf{A}^{-1}}+\|u-u_{h}\|\right)\leq\min_{\tau_{h}\in
M_{k}(\mathcal{T})}\|\sigma-\tau_{h}\|_{\mathbf{A}^{-1}}+\|u-\Pi_{k}u\|+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\sigma\|.$

(b) Suppose $\mathbf{b}_{1}=0$ and that the scalar $\gamma(x)$ is Lipschitz
continuous in $x\in\operatorname{int}(T)$, the interior of $T\in\mathcal{T}$,
with a Lipschitz constant smaller than or equal to
$\operatorname{\mathcal{L}ip}(\gamma)$. Then there exists a positive constant
$C_{2}$, which depends exclusively depends on $\beta_{0}>0$,
$\|\mathbf{A}^{1/2}\mathbf{b}_{2}\|_{L^{\infty}(\Omega)}$,
$\operatorname{\mathcal{L}ip}(\gamma)$, and the shape-regularity of
$\mathbb{T}$, such that

$C_{2}^{-1}\|\sigma-\sigma_{h}\|_{\mathbf{A}^{-1}}\leq\min_{\tau_{h}\in
M_{k}(\mathcal{T})}\|\sigma-\tau_{h}\|_{\mathbf{A}^{-1}}+\|h_{\mathcal{T}}(u-\Pi_{k}u)\|+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\sigma\|.$

The additional oscillations $\|h_{\mathcal{T}}(u-\Pi_{k}u)\|$ and
$\|h_{\mathcal{T}}(1-\Pi_{k})(\operatorname{div}\sigma)\|$ can be higher-order
contributions and then these terms explain the improved convergence of one
variant for the BDM finite element family in [8] under Assumption (A).

This article, thus, generalises earlier contributions [5], [8]-[10], [12]-[14]
for smooth or piecewise Lipschitz continuous coefficients to $L^{\infty}$
coefficients without any further assumptions. The compactness argument of
Schatz and Wang [17] for the displacement-oriented problem does not apply
immediately to the mixed formulation in $H(\operatorname{div})\times L^{2}$.
Remark $12$ below explains that no uniform $L^{2}$ approximation of the
divergence component holds. This paper therefore compensates the lack of
compactness by the computation and analysis of an optimal test function (the
dual solution $y$ in (7) of Subsection 1.4 below). Recent best-approximation
for the flux in $L^{2}$ from the medius analysis [6, 15] combines with the
compactness for the (dual) PDE. This and a careful shift of the discrete
divergence circumvents the aforementioned lack of compactness in the
divergence variable. In fact, this new methodology avoids any regularity
argument and any Fortin interpolation at all.

###  1.4 Motivation

This subsection outlines the proof of the discrete $\inf$-$\sup$ stability (6)
in an abstract framework to guide the reader through the arguments. Suppose
$X_{h}\times Y_{h}$ is a finite dimensional subspace of $H\times H$ with dual
$X_{h}^{*}\times Y_{h}^{*}$ and let $x_{h}\in S(X_{h})$, i.e., $x_{h}$ belongs
to $X_{h}$ and has norm $\|x_{h}\|_{H}=1$. Recall (5) and the well-posedness
of the problem (3). Then, the dual problem is well-posed as well and $\langle
x_{h},\bullet\rangle_{H}=b(\bullet,y)$ has a unique dual solution $y$ in the
Hilbert space $(H,\langle\bullet,\bullet\rangle_{H})$. The continuous
$\inf$-$\sup$ condition (5) shows

(7)
$\beta\,\|y\|_{H}\leq\|b(\bullet,y)\|_{H^{*}}=\|x_{h}\|_{H}=1,\quad\text{whence
}\|y\|_{H}\leq 1/\beta$

is bounded. Suppose that $y_{h}\in Y_{h}$ is a close approximation to $y$ with
$\|y-y_{h}\|_{H}\leq\varepsilon$ for some positive $\varepsilon<1/\|b\|$,
where $\|b\|$ is the operator norm of the bilinear form $b(\bullet,\bullet)$.
Since

$1=b(x_{h},y)=b(x_{h},y_{h})+b(x_{h},y-y_{h})\leq\|b(x_{h},\bullet)\|_{Y_{h}^{*}}\|y_{h}\|_{H}+\varepsilon\|b\|,$

it remains to bound $\|y_{h}\|_{H}$, e.g., with the triangle inequality

$\|y_{h}\|_{H}\leq\|y\|_{H}+\|y-y_{h}\|_{H}\leq 1/\beta+\varepsilon.$

The combination of the previous two displayed formulas gives a lower bound for
$\|b(x_{h},\bullet)\|_{Y_{h}^{*}}$. Under the assumption that $\varepsilon$ is
independent of $y$ and so of $x_{h}$, this estimate reads

(8)
$\beta\frac{1-\varepsilon\|b\|}{1+\varepsilon\beta}\leq\beta_{h}:=\inf_{x_{h}\in
X_{h}\setminus\\{0\\}}\sup_{y_{h}\in
Y_{h}\setminus\\{0\\}}\frac{b(x_{h},y_{h})}{\|x_{h}\|_{H}\|y_{h}\|_{H}}.$

This proves $\beta_{0}\leq\beta(1-\varepsilon\|b\|)/(1+\varepsilon\beta)$
provided the approximation error $\|y-y_{h}\|_{H}$ is small independently of
$X_{h}\times Y_{h}$ and $x_{h}\in S(X_{h})$. A detailed investigation in
Subsection 3.3 below reveals that the above strong form of a uniform
approximation appears neither available in the norm of
$H=H(\operatorname{div},\Omega)\times L^{2}(\Omega)$ (cf. Remark 3.5) nor
necessary for the stability under assumption (A). Recent results from a medius
analysis [6, 15] and a careful shift of the discrete divergence variable
successfully circumvent a uniform approximation in $H$.

###  1.5 Structure of the paper

Section 2 starts with the pre-compactness for uniform approximation and the
precise assumptions on the set of admissible triangulations $\mathbb{T}$. The
other two preliminary subsections concern the $L^{2}$ best-approximation of
the fluxes and some discrete approximation result for the RT finite element
family.

The stability analysis in Section 3 is based on the dual solution $y$ in the
conservative formulation characterised in Subsection 3.1. One contribution of
$y$ involves the PDE $\mathcal{L}_{2}\phi=g$ and allows for some pre-
compacness and uniform approximation in Subsection 3.2. The proof of Theorem
1.1 concludes Section 3. A combination of the stability result (6) with the
approximation arguments leads in Section 4 to Theorem 1.2, which generalises
[6, 15] to non-selfadjoint indefinite second-order linear elliptic problems.

##  2 Preliminaries

This section introduces notations used in the paper, fixes the assumptions on
the admissible triangulation $\mathbb{T}$, discusses an abstract version of
compactness argument in [17], and then recalls some $L^{2}$ best-approximation
property and concludes with an observation for the RT finite element family.

###  2.1 Notation

Standard notation on Lebesgue and Sobolev spaces $L^{2}(\Omega)$,
$L^{\infty}(\Omega)$, $H_{0}^{1}(\Omega)$, $H^{-1}(\Omega)\equiv
H_{0}^{1}(\Omega)^{*}$, and $H(\operatorname{div},\Omega)$ apply throughout
this paper. The $L^{2}$ scalar product $(\bullet,\bullet)_{L^{2}(\Omega)}$
induces the norm $\|\bullet\|:=\|\bullet\|_{L^{2}(\Omega)}$ and the
orthogonality relation $\perp$.

Whereas $\|\bullet\|$ denotes the norm in $L^{2}(\Omega)$ with the exception
of the abbreviation $\|b\|$ for the bound of the bilinear form
$b(\bullet,\bullet)$, the vector space $L^{2}(\Omega;\mathbb{R}^{n})$ is
endowed with the weighted scalar product
$(\bullet,\bullet)_{\mathbf{A}^{-1}}:=({\mathbf{A}^{-1}}\bullet,\bullet)_{L^{2}(\Omega)}$
and induced norm $\|\bullet\|_{\mathbf{A}^{-1}}:=\|\mathbf{A}^{-1/2}\bullet\|$
and so, for any $\tau\in L^{2}(\Omega;\mathbb{R}^{n})$, is its distance
$\operatorname{dist}(\tau,M_{h}):=\min_{\tau_{h}\in
M_{h}}\|\tau-\tau_{h}\|_{\mathbf{A}^{-1}}$ to any subspace $M_{h}$ of
$L^{2}(\Omega;\mathbb{R}^{n})$. The norm $\|(\tau,v)\|_{H}$ in the Hilbert
space $H(\operatorname{div},\Omega)$ is weighted with $\mathbf{A}^{-1}$ in
$L^{2}(\Omega;\mathbb{R}^{n})$ for the flux variable so the Hilbert space
$H\equiv H(\operatorname{div},\Omega)\times L^{2}(\Omega)$ has the weighted
scalar product $\langle\bullet,\bullet\rangle_{H}$ with the induced norm
$\|(\tau,v)\|_{H}$,

(9)
$\|(\tau,v)\|_{H}^{2}:=\|\tau\|_{\mathbf{A}^{-1}}^{2}+\|\operatorname{div}\tau\|^{2}+\|v\|^{2}\quad\text{for
all }(\tau,v)\in H.$

Duality brackets have the dual pairing as an index as in
$\langle\bullet,\bullet\rangle_{H^{-1}(\Omega)\times H^{1}_{0}(\Omega)}$
above. To abbreviate the definition of $\inf$-$\sup$ constants throughout this
paper, let $S(V):=\\{v\in V:\|v\|_{V}=1\\}$ for any normed linear space
$(V,\|\bullet\|_{V})$.

All emerging generic positive constants $C_{1},\dots,C_{8}$ in this paper
exclusively depend on
$\underline{\alpha},\overline{\alpha},\|\mathbf{b}\|_{L^{\infty}},\|\mathbf{b}_{1}\|_{L^{\infty}},\|\mathbf{b}_{2}\|_{L^{\infty}},$
and $\|\gamma\|_{L^{\infty}}$ as well as on $\alpha$ in (2) and $\beta$ in (5)
and on the class of admissible triangulations $\mathbb{T}$ specified in
Subsection 2.2.

###  2.2 Assumptions on the discretization

The finite element spaces are based on admissible triangulations, the set
$\mathbb{T}$ of all of those has certainly infinite cardinality; the point is
that the constants in standard interpolation error estimates become universal
through uniform shape regularity.

######  Definition 2.1 (admissible triangulations).

The set of admissible triangulations $\mathbb{T}$ is a set of shape-regular
triangulations of the polyhedral bounded Lipschitz domain
$\Omega\subset\mathbb{R}^{n}$ into simplices with uniform shape regularity and
arbitrary small mesh sizes. Let $h_{\max}(\mathcal{T}):=\max h_{\mathcal{T}}$
for the piecewise constant mesh-size $h_{\mathcal{T}}$ for
$\mathcal{T}\in\mathbb{T}$, defined by $h_{\mathcal{T}}|_{T}:=\text{\rm
diam}(T)$ in $T\in\mathcal{T}$, and abbreviate
$\mathcal{T}(\delta):=\\{\mathcal{T}\in\mathbb{T}:h_{\max}(\mathcal{T})\leq\delta\\}$.

Given $\mathcal{T}\in\mathbb{T}$, let $P_{k}(T)$ denote the polynomials of
total degree at most $k\in\mathbb{N}_{0}$ seen as functions on
$T\in\mathcal{T}\in\mathbb{T}$ and set $P_{k}(\mathcal{T}):=\\{v_{k}\in
L^{\infty}(\Omega):\forall T\in\mathcal{T},\;v_{k}|_{T}\in P_{k}(T)\\}$. Let
$\Pi_{k}:L^{2}(\Omega)\to L^{2}(\Omega)$ be the $L^{2}$ projection onto
$P_{k}(\mathcal{T})$ with respect to $\mathcal{T}\in\mathbb{T}$.

######  Definition 2.2 (discrete spaces).

Any $\mathcal{T}\in\mathbb{T}$ is associated to the finite-dimensional
subspace $V(\mathcal{T})=M_{k}(\mathcal{T})\times P_{k}(\mathcal{T})$ of
$V:=L^{2}(\Omega;\mathbb{R}^{n})\times L^{2}(\Omega)$ with
$M_{k}(\mathcal{T}):=RT_{k}(\mathcal{T})$ or
$M_{k}(\mathcal{T}):=BDM_{k}(\mathcal{T})$ of order $k\in\mathbb{N}_{0}$ from
[2].

The best-approximation error reads
$\operatorname{dist}(v,V(\mathcal{T})):=\inf\\{\|v-v_{h}\|:v_{h}\in
V(\mathcal{T})\\}$ with the weighted $L^{2}$ norm,
$\|v\|^{2}=\|\tau\|^{2}_{\mathbf{A}^{-1}}+\|w\|^{2}$ for $v=(\tau,w)\in V$.
The density of smooth functions and standard approximation results for smooth
functions proves the well-known pointwise convergence in the sense that each
$v\in V$ satisfies [2]

(10) $\lim_{\delta\to
0^{+}}\sup_{\mathcal{T}\in\mathbb{T}(\delta)}\operatorname{dist}(v,V(\mathcal{T}))=0.$

###  2.3 Pre-compactness

This subsection adopts the key argument of [17].

######  Lemma 2.3 (uniform approximation on compact sets).

Suppose that $K$ is a non-empty pre-compact subset of
$V:=L^{2}(\Omega;\mathbb{R}^{n})\times L^{2}(\Omega)$ with (10) for each $v\in
K$. Then

(11) $\lim_{\delta\to 0^{+}}\sup_{v\in
K}\sup_{\mathcal{T}\in\mathbb{T}(\delta)}\operatorname{dist}(v,V(\mathcal{T}))=0.$

######  Proof 2.4.

Given any $\varepsilon>0$ and $v\in K$, let $B(v,\varepsilon/2)$ be the open
ball in $V$ with center $v$ and radius $\varepsilon/2$. The open cover
$\\{B(v,\varepsilon/2):v\in K\\}$ of the compact set $\overline{K}$ contains a
finite sub-cover and so there exist $k_{1},\dots,k_{J}\in K$ with
$K\subset\bigcup_{j=1,\dots,J}B(k_{j},\varepsilon/2)$. For each $k_{j}\in K$,
(10) leads to $\delta_{j}>0$ such that $\mathcal{T}\in\mathcal{T}(\delta_{j})$
implies $\operatorname{dist}(k_{j},V(\mathcal{T}))<\varepsilon/2$. Then
$\delta:=\min\\{\delta_{1},\dots,\delta_{J}\\}$ implies
$\mathbb{T}(\delta)\subset\bigcap_{j=1,\dots,J}\mathbb{T}(\delta_{j})$. Given
any $\mathcal{T}\in\mathbb{T}(\delta)$ and any $v\in
K\subset\bigcup_{j=1,\dots,J}B(k_{j},\varepsilon/2)$, there exists
$j\in\\{1,\dots,J\\}$ with $\|v-k_{j}\|<\varepsilon/2$. Since
$\mathcal{T}\in\mathbb{T}(\delta_{j})$,
$\operatorname{dist}(k_{j},V(\mathcal{T}))<\varepsilon/2$. This and a triangle
inequality show
$\operatorname{dist}(v,V(\mathcal{T}))\leq\|v-k_{j}\|+\operatorname{dist}(k_{j},V(\mathcal{T}))<\varepsilon/2+\varepsilon/2=\varepsilon.$

The application of the previous lemma to the finite element approximation of
the solution of the PDE reads as follows.

######  Lemma 2.5 (uniform approximation of solutions).

For any $\varepsilon>0$ there exists some $\delta>0$ such that, given any
$g\in L^{2}(\Omega)$ and the weak solution $\phi\in H^{1}_{0}(\Omega)$ to
$\mathcal{L}_{2}\phi=g$ (with $\mathcal{L}_{2}$ from (1)), the vector
$v:=(\mathbf{A}\nabla\phi,\mathbf{b}\cdot\nabla\phi+\gamma\phi)\in V$
satisfies
$\sup_{\mathcal{T}\in\mathbb{T}(\delta)}\operatorname{dist}(v,V(\mathcal{T}))\leq\epsilon\,\|g\|.$

######  Proof 2.6.

The linear and bounded bijective differential operator
$\mathcal{L}_{2}:H^{1}_{0}(\Omega)\to H^{-1}(\Omega)$ has a bounded inverse.
The embedding $\iota:L^{2}(\Omega)\hookrightarrow H^{-1}(\Omega)$ is compact
and so is the composition $\mathcal{L}_{2}^{-1}\circ\iota:L^{2}(\Omega)\to
H^{1}_{0}(\Omega)$. Define the operator $T:L^{2}(\Omega)\to
L^{2}(\Omega;\mathbb{R}^{n})\times L^{2}(\Omega)$ for any $g\in L^{2}(\Omega)$
by

$T(g):=(\mathbf{A}\nabla\phi,\mathbf{b}\cdot\nabla\phi+\gamma\phi)\quad\text{with}\quad\phi:=\mathcal{L}_{2}^{-1}g.$

Since $\mathcal{L}_{2}^{-1}\circ\iota$ is compact, $K:=T(S(L^{2}(\Omega))$ is
pre-compact in $V=L^{2}(\Omega;\mathbb{R}^{n})\times L^{2}(\Omega)$. Given any
$\varepsilon>0$ the approximation result (10) and Lemma 2.3 lead to a positive
$\delta$ with (11). Consequently, the assertion
$\sup_{\mathcal{T}\in\mathbb{T}(\delta)}\operatorname{dist}(v,V(\mathcal{T}))\leq\epsilon\,\|g\|$
holds for all $g\in S(L^{2}(\Omega))$ and corresponding
$v:=(\mathbf{A}\nabla\phi,\mathbf{b}\cdot\nabla\phi+\gamma\phi)\in V$. A
rescaling proves the result for all $g\in L^{2}(\Omega)$.

###  2.4 $L^{2}$ best-approximation of the fluxes

The medius analysis of mixed finite element methods employs arguments from a
priori and a posteriori error analysis [6, 15] to prove new $L^{2}$ best-
approximation results. Recall that $\Pi_{k}$ is the $L^{2}$ projection onto
$P_{k}(\mathcal{T})$ and $h_{\mathcal{T}}$ is the mesh-size associated to
$\mathcal{T}$.

######  Lemma 2.7 (flux $L^{2}$ best-approximation).

There exists a constant $C_{3}$, which depends on the shape-regularity in
$\mathcal{T}$, on $\Omega$ and on $\underline{\alpha},\overline{\alpha}$, such
for any $\mathbf{p}\in H(\operatorname{div},\Omega)$ and any
$\mathcal{T}\in\mathbb{T}$, there exists $\mathbf{p}_{h}\in
M_{k}(\mathcal{T})$ such that
$\operatorname{div}\mathbf{p}_{h}=\Pi_{k}\operatorname{div}\mathbf{p}$ and

$C_{3}^{-1}\|\mathbf{p}-\mathbf{p}_{h}\|_{\mathbf{A}^{-1}}\leq\operatorname{dist}(\mathbf{p},M_{k}(\mathcal{T}))+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\mathbf{p}\|.$

This is the $L^{2}$ best-approximation result from [15, Lemma 5.1] for mixed
finite element approximations for the unit matrix $\mathbf{A}$. Although with
a different focus, the paper [6] introduces a general framework with a mesh-
dependent norm $\|\bullet\|_{h}$ in $P_{k}(\mathcal{T})$; while [11, Eq (3.6)]
presents a localized refinement of this lemma.

######  Proof 2.8.

Given $\mathbf{p}\in H(\operatorname{div},\Omega)$, the right-hand sides
$F(w):=(w,\operatorname{div}\mathbf{p})_{L^{2}(\Omega)}$ and
$G(\mathbf{q}):=(\mathbf{p},\mathbf{q})$ lead in the elliptic mixed
formulation (for the Laplacian)

$(\sigma,\mathbf{q})-(u,\operatorname{div}\mathbf{q})_{L^{2}(\Omega)}+(w,\operatorname{div}\mathbf{p})_{L^{2}(\Omega)}=G(\mathbf{q})+F(w)\quad\text{for
all }(\mathbf{q},w)\in H$

to the unique solution $(\sigma,u)\equiv(\mathbf{p},0)\in H$. Its straight-
forward mixed finite element discretisation substitutes $H$ by
$V_{h}:=M_{k}(\mathcal{T})\times P_{k}(\mathcal{T})$ and leads to a unique
discrete solution $(\mathbf{p}_{h},v_{h})\in V_{h}$ with
$\operatorname{div}\mathbf{p}_{h}=\Pi_{k}\operatorname{div}\mathbf{p}$. This
and [6, Thm 2.2] lead to the asserted best-approximation result (in terms of
(non-weighted) $L^{2}$ norms)

$C_{4}^{-1}\|\mathbf{p}-\mathbf{p}_{h}\|\leq\inf_{{\mathbf{q}}_{h}\in
M_{k}(\mathcal{T})}\|\mathbf{p}-{\mathbf{q}}_{h}\|+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\mathbf{p}\|.$

The constant $C_{4}$ from [6, 15] does not depend on the coefficients ${\bf
A},$ b, $\gamma$ but depends on the shape-regularity in $\mathcal{T}$ and on
$\Omega$. The equivalence of norms concludes the proof and leads to the
asserted constant $C_{3}$, which depends on $C_{4}$ and
$\underline{\alpha},\overline{\alpha}$.

###  2.5 A discrete approximation result for Raviart-Thomas functions

In any space-dimension $n$ and degree $k$, the RT functions satisfy a rather
particular approximation estimate with the componentwise $L^{2}$ projection
$\Pi_{k}$ onto $P_{k}(\mathcal{T})$.

######  Lemma 2.9.

Any $\tau_{RT}\in RT_{k}(\mathcal{T})$ satisfies
$\|\tau_{RT}-\Pi_{k}\tau_{RT}\|\leq\frac{n}{(n+1)(n+k)}\|h_{\mathcal{T}}\operatorname{div}\tau_{RT}\|.$

The proof will be postponed to the appendix because of its focus on the RT
finite element shape functions. The statement of the above lemma fails for the
BDM finite element family.

##  3 Stability analysis

This section deals with approximation of fluxes and stability result. The
design of a test function in the proof of a discrete $\inf$-$\sup$ condition
is based on the characterisation and approximation of a dual solution.

###  3.1 Dual solution and conservative formulation

The inner structure of the dual solution $y$ exploits the elliptic PDE and
generates some compactness argument in the subsequent subsection. Recall that
the operator $\mathcal{L}_{2}:H^{1}_{0}(\Omega)\to H^{-1}(\Omega)$ from (1) is
bijective.

######  Theorem 3.1 (dual solution in conservative formulation).

Suppose $\mathbf{b}_{1}:=\mathbf{A}^{-1}\mathbf{b}$ and $\mathbf{b}_{2}\equiv
0$ a.e. in $\Omega$ in (4). Then $x=(\sigma,u)\in H$ and $y=(\zeta,z)\in H$
satisfy $\langle x,\bullet\rangle_{H}=b(\bullet,y)$ in $H$ if and only if

$\zeta=\sigma-\mathbf{A}\nabla\phi\quad\text{and}\quad
z=\operatorname{div}\sigma-\phi\quad\text{ a.e. in }\Omega$

for the weak solution $\phi\in H^{1}_{0}(\Omega)$ to
$\mathcal{L}_{2}\phi=g:=\mathbf{b}\cdot\mathbf{A}^{-1}\sigma+(\gamma-1)\,\operatorname{div}\sigma-u\in
L^{2}(\Omega)$.

The function $\phi$ originates from a known $L^{2}$ orthogonal decomposition

(12) $L^{2}(\Omega;\mathbb{R}^{n})=\nabla H^{1}_{0}(\Omega)\oplus H(\rm
div{=}0,\Omega)$

with $H(\rm div{=}0,\Omega):=\\{\tau\in
H(\operatorname{div},\Omega):\,\operatorname{div}\tau=0\text{ a.e. in
}\Omega\\}$. The decomposition (12) is also useful in the proof of equivalence
of the displacement formulation with the differential operators in (1) to the
mixed formulations with (3).

Proof of Theorem 3.1. For the general version of the bilinear form
$b(\bullet,\bullet)$, the equation $\langle x,\bullet\rangle_{H}=b(\bullet,y)$
is equivalent to
$u=\mathbf{b}_{1}\cdot\zeta+\gamma\,z-\operatorname{div}\zeta$ a.e. in
$\Omega$ and

(13)
$\displaystyle(\tau,\mathbf{A}^{-1}(\zeta-\sigma)-z\mathbf{b}_{2})_{L^{2}(\Omega)}+(z-\operatorname{div}\sigma,\operatorname{div}\tau)_{L^{2}(\Omega)}=0\quad\text{for
all }\tau\in H(\operatorname{div},\Omega).$

The test with $\tau\in H(\rm div{=}0,\Omega)$ proves that
$\mathbf{A}^{-1}(\sigma-\zeta)+z\mathbf{b}_{2}\perp H(\rm div{=}0,\Omega)$ and
so (12) leads to $\phi\in H^{1}_{0}(\Omega)$ with

$\mathbf{A}\nabla\phi=\sigma-\zeta+z\mathbf{A}\mathbf{b}_{2}\quad\text{ a.e.
in }\Omega.$

This identity allows the substitution of
$\mathbf{A}^{-1}(\zeta-\sigma)-z\mathbf{b}_{2}$ in the above formula (13) with
general $\tau\in H(\operatorname{div},\Omega)$. Then, an integration by parts
shows the resulting identity
$(\phi+z-\operatorname{div}\sigma,\operatorname{div}\tau)=0$. The surjectivity
of $\operatorname{div}:H(\operatorname{div},\Omega)\to L^{2}(\Omega)$ proves

$\operatorname{div}\sigma=\phi+z\quad\text{ a.e. in }\Omega.$

The combination of the three preceding identities leads to the PDE

$-\operatorname{div}(\mathbf{A}\nabla\phi)+\mathbf{b}_{1}\cdot\mathbf{A}\nabla\phi+\gamma\phi=-\operatorname{div}(z\mathbf{A}\mathbf{b}_{2})+(\mathbf{b}_{1}\cdot\mathbf{A}\mathbf{b}_{2})\,z+\mathbf{b}_{1}\cdot\sigma+(\gamma-1)\,\operatorname{div}\sigma-u$

in the sense of distributions. Since $\mathbf{b}_{2}=0$, the right-hand side
$g$ belongs to $L^{2}$. This proves one direction of the assertion; the direct
proof of the converse is omitted. $\Box$

######  Remark 3.2 (no divergence formulation).

The proof shows the extra term
$-\operatorname{div}(z\mathbf{A}\mathbf{b}_{2})\in H^{-1}(\Omega)$ in case (4)
is considered for non-zero $\mathbf{b}_{2}\in
L^{\infty}(\Omega;\mathbb{R}^{n})$. This term does not belong to
$L^{2}(\Omega)$ under Assumption (A) and is, therefore, excluded.

###  3.2 Approximation of the fluxes

The subsequent lemma describes the uniform approximation of the flux variable
by a combination of the compactness argument and the $L^{2}$ best-
approximation of Subsections 2.3 and 2.4.

######  Lemma 3.3 (flux approximation).

Given any $\varepsilon>0$, there exists $\delta>0$ such that the following
holds for all $\mathcal{T}\in\mathbb{T}(\delta)$ and $g\in L^{2}(\Omega)$.
There exists some $\mathbf{p}_{h}\in M_{k}(\mathcal{T})$ that approximates
$\mathbf{p}:=\mathbf{A}\nabla\phi\in H(\operatorname{div},\Omega)$ for the
weak solution $\phi\in H^{1}_{0}(\Omega)$ to $\mathcal{L}_{2}\phi=g$ with

$\operatorname{div}\mathbf{p}_{h}=\Pi_{k}\operatorname{div}\mathbf{p}\quad\text{and}\quad\|\mathbf{p}-\mathbf{p}_{h}\|_{\mathbf{A}^{-1}}\leq\epsilon\|g\|.$

######  Proof 3.4.

Given any $\varepsilon>0$ and the constant $C_{3}$ from Lemma 2.7, Lemma 2.5
leads to a positive $\delta\leq\min\\{1,2^{-1}\epsilon/C_{3}\\}$ with

$\sup_{\mathcal{T}\in\mathbb{T}(\delta)}\operatorname{dist}\left((\mathbf{A}\nabla\phi,\mathbf{b}\cdot\nabla\phi+\gamma\phi),V(\mathcal{T})\right)\leq
2^{-3/2}\epsilon/C_{3}\,\|g\|$

(the distance is with respect to the weighted norm
$\|\bullet\|_{\mathbf{A}^{-1}}$ in $L^{2}(\Omega;\mathbb{R}^{n})$ and
$\|\bullet\|$ in $L^{2}(\Omega)$). Lemma 2.7 applies to
$\mathbf{p}:=\mathbf{A}\nabla\phi$ with
$\operatorname{div}\mathbf{p}=\mathbf{b}\cdot\nabla\phi+\gamma\phi-g\in
L^{2}(\Omega)$ and, for any $\mathcal{T}\in\mathbb{T}(\delta)$, leads to some
approximation $\mathbf{p}_{h}\in M_{k}(\mathcal{T})$ with
$\operatorname{div}\mathbf{p}_{h}=\Pi_{k}\operatorname{div}\mathbf{p}$ and

$\displaystyle C_{3}^{-1}\|\mathbf{p}-\mathbf{p}_{h}\|_{\mathbf{A}^{-1}}$
$\displaystyle\leq\operatorname{dist}(\mathbf{p},M_{k}(\mathcal{T}))+\delta\|(1-\Pi_{k})(\mathbf{b}\cdot\nabla\phi+\gamma\phi-g)\|$
$\displaystyle\leq\operatorname{dist}(\mathbf{p},M_{k}(\mathcal{T}))+\operatorname{dist}(\mathbf{b}\cdot\nabla\phi+\gamma\phi,P_{k}(\mathcal{T}))+\delta\,\|g\|$
$\displaystyle\leq
2^{1/2}\operatorname{dist}((\mathbf{p},\mathbf{b}\cdot\nabla\phi+\gamma\phi),V(\mathcal{T}))+\delta\,\|g\|\leq
C_{3}^{-1}\epsilon\,\|g\|.$

This concludes the proof.

######  Remark 3.5 (no uniform approximation in $H(\operatorname{div})$).

Lemma 3.3 does not state a uniform approximation estimate for the divergence
and, in fact, an estimate of the form
$\|\operatorname{div}(\mathbf{p}-\mathbf{p}_{h})\|\leq\epsilon\|g\|$ cannot
hold in general. To see this, adopt the notation of the proof of Lemma 3.3 and
a reverse triangle inequality for

$\|g-\Pi_{k}g\|-\|\operatorname{div}(\mathbf{p}-\mathbf{p}_{h})\|\leq\|(1-\Pi_{k})(\mathbf{b}\cdot\nabla\phi+\gamma\phi)\|\leq\epsilon/(2C_{3})\,\|g\|.$

The flux $\mathbf{p}=\mathbf{A}\nabla\mathcal{L}_{2}^{-1}g$ depends on $g\in
S(L^{2}(\Omega))$ and so does the crucial term
$\|\operatorname{div}(\mathbf{p}-\mathbf{p}_{h})\|=\|(1-\Pi_{k})\operatorname{div}\mathbf{p}\|$.
Hence, $\sup\\{\|g-\Pi_{k}g\|:g\in S(L^{2}(\Omega))\\}=1$ implies

$1-\epsilon/(2C_{3})\leq\sup\\{\|(1-\Pi_{k})\operatorname{div}(\mathbf{A}\nabla\mathcal{L}_{2}^{-1}g)\|:g\in
S(L^{2}(\Omega))\\}.$

Therefore, the approximation error
$\|\operatorname{div}(\mathbf{p}-\mathbf{p}_{h})\|$ will not tend to zero
uniformly for all $g\in S(L^{2}(\Omega))$ as $\epsilon$ and $\delta$ tend to
zero. $\Box$

######  Example 3.6 (RT approximation in $H(\operatorname{div})$ for
particular $g$).

The stability analysis in Subsection 1.4 concerns a discrete
$x_{h}:=(\sigma_{h},u_{h})$ with norm $\|(\sigma_{h},u_{h})\|_{H}=1$ and leads
in Theorem 3.1 to the particular right-hand side
$g:=\mathbf{b}\cdot\mathbf{A}^{-1}\sigma_{h}+(\gamma-1)\,\operatorname{div}\sigma_{h}-u_{h}$
with $\|g\|\leq C_{5}<\infty$ for the essential supremum $C_{5}^{2}$ of
$|\mathbf{A}^{-1/2}\mathbf{b}|^{2}+|\gamma-1|^{2}+1$ in $\Omega$. This $g$
allows for a uniform approximation of $\mathbf{p}$ by $\mathbf{p}_{h}$ in
$H(\operatorname{div},\Omega)$ for the RT finite element family.

For instance, in the extreme case of piecewise constant coefficients,
$g-\Pi_{k}g=\mathbf{b}\cdot\mathbf{A}^{-1}(1-\Pi_{k})\sigma_{h}$. With
$C_{6}:=\|\mathbf{A}^{-1}\mathbf{b}\|_{L^{\infty}(\Omega)}$, Lemma 2.9 shows
$\|g-\Pi_{k}g\|\leq\delta C_{6}$. The combination with Lemma 3.3 lead to
$\mathbf{p}_{h}$ with
$\|\mathbf{p}-\mathbf{p}_{h}\|_{H(\operatorname{div},\Omega)}\leq\epsilon\,C_{5}+\delta\,C_{6}.$
This and the arguments of Subsection 1.4 lead to the discrete stability (6).

###  3.3 Proof of Theorem 1.1

Given any $0<\varepsilon<\beta/\|b\|$, choose $\delta>0$ as in Lemma 3.3.
Suppose $\mathcal{T}\in\mathbb{T}(\delta)$ and let
$x_{h}=(\sigma_{h},u_{h})\in V_{h}:=M_{k}(\mathcal{T})\times
P_{k}(\mathcal{T})$ have norm $\|x_{h}\|_{H}=1$ and define
$g:=\mathbf{b}\cdot\mathbf{A}^{-1}\sigma_{h}+(\gamma-1)\,\operatorname{div}\sigma_{h}-u_{h}$.
Replace $x$ by $x_{h}$ in Theorem 3.1 and let $\phi\in H^{1}_{0}(\Omega)$
solve $\mathcal{L}_{2}\phi=g$ to define
$\zeta=\sigma_{h}-\mathbf{A}\nabla\phi$ and
$z=\operatorname{div}\sigma_{h}-\phi$. Then $y=(\zeta,z)\in H$ is the dual
solution and solves $\langle x_{h},\bullet\rangle_{H}=b(\bullet,y)$ in $H$ for
the bilinear form (4) (with $\mathbf{b}_{1}:=\mathbf{A}^{-1}\mathbf{b}$ and
$\mathbf{b}_{2}\equiv 0$ a.e. in $\Omega$). Lemma 3.3 applies to
$\mathbf{p}:=\mathbf{A}\nabla\phi$ and leads to $\mathbf{p}_{h}$ with
$\operatorname{div}\mathbf{p}_{h}=\Pi_{k}\operatorname{div}\mathbf{p}$ and
$\|\mathbf{p}-\mathbf{p}_{h}\|_{\mathbf{A}^{-1}}\leq\varepsilon\,C_{5}$.
Hence, $\zeta_{h}:=\sigma_{h}-\mathbf{p}_{h}$ and $z_{h}:=\Pi_{k}z$ define
$y_{h}:=(\zeta_{h},z_{h})\in V_{h}$ with

(14)
$\|y_{h}\|_{H}\leq\|y\|_{H}+\|\zeta-\zeta_{h}\|_{\mathbf{A}^{-1}}\leq\|y\|_{H}+\varepsilon\,C_{5}$

as
$\|\operatorname{div}\zeta_{h}\|=\|\Pi_{k}\operatorname{div}\zeta\|\leq\|\operatorname{div}\zeta\|$,
$\|z_{h}\|=\|\Pi_{k}z\|\leq\|z\|$ and $\operatorname{div}\sigma_{h}\in
P_{k}(\mathcal{T})$. Moreover,

$\|\zeta-\zeta_{h}\|_{\mathbf{A}^{-1}}^{2}+\|z-z_{h}\|^{2}\leq\varepsilon^{2}\,C_{5}^{2}+\|\phi-\Pi_{k}\phi\|^{2}.$

Piecewise Poincaré inequalities (with the Payne-Weinberger constant $1/\pi$
for convex domains [16]) show
$\|\phi-\Pi_{k}\phi\|\leq\delta/\pi\|\nabla\phi\|$. Recall that
$H^{1}_{0}(\Omega)$ is endowed with the seminorm $\|\nabla\bullet\|$ and let
$C_{F}$ denote the constant in the Friedrichs inequality $\|\bullet\|\leq
C_{F}\|\nabla\bullet\|$ in $H^{1}_{0}(\Omega)$. Note that (2) leads to
$\alpha\,\|\nabla\phi\|\leq\sup\\{a(\psi,\phi):\psi\in
H^{1}_{0}(\Omega),\;\|\nabla\psi\|=1\\}$. Hence
$a(\psi,\phi)=\langle\mathcal{L}_{2}\phi,\psi\rangle_{H^{-1}(\Omega)\times
H^{1}_{0}(\Omega)}=\langle g,\psi\rangle_{H^{-1}(\Omega)\times
H^{1}_{0}(\Omega)}\leq C_{F}\|g\|$ implies $\alpha\,\|\nabla\phi\|\leq
C_{F}\,C_{5}$. The combination with Poincaré inequalities shows
$\|\phi-\Pi_{k}\phi\|\leq\delta C_{F}\,C_{5}/(\alpha\,\pi)$ and so

$\|\zeta-\zeta_{h}\|_{\mathbf{A}^{-1}}^{2}+\|z-z_{h}\|^{2}\leq\varepsilon^{2}\,C_{5}^{2}+\delta^{2}C_{F}^{2}\,C_{5}^{2}/(\alpha^{2}\,\pi^{2})=:(\varepsilon^{\prime})^{2}.$

Since $u_{h}\perp\operatorname{div}(\zeta-\zeta_{h})$ and
$z-z_{h}\perp\operatorname{div}\sigma_{h}$ ($\perp$ denotes orthogonality in
$L^{2}(\Omega)$),

$\displaystyle b(x_{h},y-y_{h})$
$\displaystyle=(\sigma_{h}+\mathbf{b}u_{h},\zeta-\zeta_{h})_{\mathbf{A}^{-1}}+(\gamma\,u_{h},z-z_{h})$
$\displaystyle\leq\varepsilon^{\prime}\left(\|\sigma_{h}+\mathbf{b}u_{h}\|^{2}_{\mathbf{A}^{-1}}+\|(1-\Pi_{k})(\gamma\,u_{h})\|^{2}\right)^{1/2}$
$\displaystyle\leq\varepsilon^{\prime}\left(2\|\sigma_{h}\|^{2}_{\mathbf{A}^{-1}}+2\|\mathbf{A}^{-1/2}\mathbf{b}\|_{L^{\infty}(\Omega)}^{2}\|u_{h}\|^{2}+\|\gamma\|_{L^{\infty}(\Omega)}^{2}\|u_{h}\|^{2}\right)^{1/2}\leq\varepsilon^{\prime}\,C_{7}$

with the constant
$C_{7}^{2}:=\max\\{2,2\|\mathbf{A}^{-1/2}\mathbf{b}\|_{L^{\infty}(\Omega)}^{2}+\|\gamma\|_{L^{\infty}(\Omega)}^{2}\\}$.
The arguments of Subsection 1.4 lead to

(15)
$1=\|x_{h}\|^{2}_{H}=b(x_{h},y)=b(x_{h},y_{h})+b(x_{h},y-y_{h})\leq\|b(x_{h},\bullet)\|_{V_{h}^{*}}\|y_{h}\|_{H}+\varepsilon^{\prime}\,C_{7}.$

Since $\beta\|y\|_{H}\leq\|b(\bullet,y)\|_{H^{*}}=\|x_{h}\|_{H}=1$ implies
$\|y\|_{H}\leq 1/\beta$, (14) reads $\|y_{h}\|_{H}\leq
1/\beta+\varepsilon\,C_{5}$. This and (15) verify

$1-\varepsilon^{\prime}\,C_{7}\leq\|b(x_{h},\bullet)\|_{V_{h}^{*}}\|y_{h}\|_{H}\leq\|b(x_{h},\bullet)\|_{V_{h}^{*}}(1/\beta+\varepsilon\,C_{5}).$

Since $x_{h}$ was arbitrary in $S(V_{h})$ (with $V_{h}$ endowed with the norm
in $H$),

$\beta\frac{1-\varepsilon^{\prime}\,C_{7}}{1+\varepsilon\,\beta\,C_{5}}\leq\beta_{h}:=\inf_{x_{h}\in
S(V_{h})}\sup_{y_{h}\in S(V_{h})}b(x_{h},y_{h}).$

Relabelling $\epsilon$ and $\delta$ proves the assertion: For any
$0<\beta_{0}<\beta$, there exists $\delta>0$ with (6). This establishes the
theorem for the conservative version of the mixed finite element
discretisation with $\mathbf{b}_{1}:=\mathbf{A}^{-1}\mathbf{b}$ and
$\mathbf{b}_{2}\equiv 0$ a.e. in $\Omega$.

To deduce the same $\inf$-$\sup$ constant for the other variant, define
$b_{1}(\bullet,\bullet)$ (resp. $b_{2}(\bullet,\bullet)$) by (6) for
$\mathbf{b}_{1}:=\mathbf{A}^{-1}\mathbf{b}$ and $\mathbf{b}_{2}\equiv 0$
(resp. ${\mathbf{b}}_{1}:=0$ and
${\mathbf{b}}_{2}:={\mathbf{b}}\cdot\mathbf{A}^{-1}$) a.e. in $\Omega$. The
above proof shows $0<\beta_{0}\leq\beta_{h}$ and it is elementary to see that

$\displaystyle\beta_{h}$ $\displaystyle=\inf_{(\tau_{h},v_{h})\in
S(V_{h})}\sup_{(\sigma_{h},u_{h})\in
S(V_{h})}b_{1}((\tau_{h},-v_{h}),(\sigma_{h},-u_{h})).$

A direct calculation shows
$b_{1}((\tau,-v),(\sigma,-u))=b_{2}((\sigma,u),(\tau,v))$ for all
$(\sigma,u),(\tau,v)\in H$. This and a duality argument (singular values of a
square matrix coincide with those of its transposed) in the last equality show

$\beta_{h}=\inf_{(\tau_{h},v_{h})\in S(V_{h})}\sup_{(\sigma_{h},u_{h})\in
S(V_{h})}b_{2}((\sigma_{h},u_{h}),(\tau_{h},v_{h}))=\inf_{x_{h}\in
S(V_{h})}\sup_{y_{h}\in S(V_{h})}b_{2}(x_{h},y_{h}).$

Hence, the divergence formulation has the same discrete $\inf$-$\sup$ constant
$\beta_{h}$. $\Box$

######  Remark 3.7 ($\delta$ dependence).

The size of $\delta$ in (6) is hidden behind a compactness argument of Lemma
3.3. Besides the norms and parameters mentioned in Assumption (A), the mapping
properties of $\mathcal{L}_{2}^{-1}$ are of relevance as well. A review of the
proofs of this paper shows that there is a finite sub-cover of
$S(L^{2}(\Omega))$ with small balls in $H^{-1}(\Omega)$ that leads to a finite
number of (without loss of generality) smooth functions $k_{1},\dots,k_{J}$ as
in the proof of Lemma 2.3. The size of $\delta$ is related to the
approximation properties of the weak solutions $\Phi_{j}$ to
$\mathcal{L}_{2}\Phi_{j}=k_{j}$ a.e. The regularity properties of $\Phi_{j}\in
H^{1}_{0}(\Omega)$ are not characterised for Assumption (A): In fact, it is
unknown whether $\Phi_{j}$ belongs to any $H^{1+s}(\Omega)$ for any $s>0$.
Under Assumption (B) and reduced elliptic regularity, however, the afore
mentioned approximation properties could be quantified more and reveal further
information on $\delta$.

##  4 $L^{2}$ Best-approximation

The notation of Theorem 1.2 applies throughout this section with continuous
and discrete solutions $x=(\sigma,u)$ and $x_{h}=(\sigma_{h},u_{h})$.

###  4.1 Proof of Theorem 1.2.a

Given $\mathbf{p}:=\sigma\in H(\operatorname{div},\Omega)$ and
$\mathcal{T}\in\mathbb{T}(\delta)$, choose $\sigma_{h}^{*}:=\mathbf{p}_{h}\in
M_{k}(\mathcal{T})$ as in Lemma 2.7, and define
$e_{h}:=(\sigma_{h}-\sigma_{h}^{*},u_{h}-\Pi_{k}u)\in V_{h}$. Given
$\beta_{0}>0$ in (6) there exists some $y_{h}=(\tau_{h},v_{h})\in V_{h}$ with
$\|y_{h}\|_{H}=1$ and

$\beta_{0}\,\|e_{h}\|_{H}\leq
b(e_{h},y_{h})=b((\sigma_{h}-\sigma_{h}^{*},u_{h}-\Pi_{k}u),y_{h})=b((\sigma-\sigma_{h}^{*},u-\Pi_{k}u),y_{h}).$

Since $u-\Pi_{k}u\perp\operatorname{div}\tau_{h}$ and
$v_{h}\perp\operatorname{div}(\sigma-\sigma^{*}_{h})$, the last term is equal
to

$(\sigma-\sigma^{*}_{h},\tau_{h}-v_{h}\,\mathbf{A}\mathbf{b}_{2})_{\mathbf{A}^{-1}}+(u-\Pi_{k}u,\mathbf{b}_{1}\cdot\tau_{h}+\gamma\,v_{h})_{L^{2}(\Omega)}\leq
C_{8}\,\|(\sigma-\sigma^{*}_{h},u-\Pi_{k}u)\|_{L}$

in terms of the weighted $L^{2}$ norm $\|\bullet\|_{L}\leq\|\bullet\|_{H}$ in
$H$ with
$\|(\tau_{h},v_{h})\|_{L}^{2}:=\|\tau_{h}\|_{\mathbf{A}^{-1}}^{2}+\|v_{h}\|^{2}\leq
1$ and with
$C_{8}^{2}=1+\|\mathbf{A}^{1/2}\mathbf{b}_{1}\|_{L^{\infty}(\Omega)}^{2}+\|\mathbf{A}^{1/2}\mathbf{b}_{2}\|_{L^{\infty}(\Omega)}^{2}+\|\gamma\|_{L^{\infty}(\Omega)}^{2}$.
Consequently,
$\|e_{h}\|_{H}\leq\Big{(}C_{8}/\beta_{0}\Big{)}\,\|(\sigma-\sigma^{*}_{h},u-\Pi_{k}u)\|_{L}$.
Lemma 2.7 shows

(16)
$C_{3}^{-1}\|\sigma-\sigma^{*}_{h}\|_{\mathbf{A}^{-1}}\leq\operatorname{dist}(\sigma,M_{k}(\mathcal{T}))+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\sigma\|.$

This and the distance $\operatorname{dist}_{L}$ (measured in the norm
$\|\bullet\|_{L}$) lead to

$\|e_{h}\|_{H}\leq
C_{8}\max\\{1,C_{3}\\}/\beta_{0}\,\left(\operatorname{dist}_{L}((\sigma,u),V_{h})+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\sigma\|\right).$

This, (16), and a triangle inequality conclude the proof. $\Box$

###  4.2 Proof of Theorem 1.2.b

Throughout this subsection, let $\mathbf{b}_{1}\equiv 0$ and
$\mathbf{b}_{2}:=\mathbf{A}^{-1}\mathbf{b}$ a.e. in $\Omega$ in (4) and let
$f\in L^{2}(\Omega)$ be a fixed right-hand side for the continuous and
discrete problem $\mathcal{L}_{2}u=f$.

Return to the proof of the previous subsection with $e_{h}$ and follow the
first lines until

$\beta_{0}\,\|e_{h}\|_{H}\leq
b(e_{h},y_{h})=(\sigma-\sigma^{*}_{h},\tau_{h}-v_{h}\,\mathbf{A}\mathbf{b}_{2})_{\mathbf{A}^{-1}}+(u-\Pi_{k}u,\gamma\,v_{h})_{L^{2}(\Omega)}.$

Abbreviate $\overline{\gamma}:=\Pi_{0}\gamma$ and utilize the Lipschitz
continuity of the coefficients $\gamma$ on each simplex

$\|\Pi_{k}(\gamma(u-\Pi_{k}u))\|=\|\Pi_{k}((\gamma-\overline{\gamma})(u-\Pi_{k}u))\|\leq\operatorname{\mathcal{L}ip}(\gamma)\|h_{\mathcal{T}}(u-\Pi_{k}u)\|.$

This controls the above term
$(u-\Pi_{k}u,\gamma\,v_{h})_{L^{2}(\Omega)}\leq\|\Pi_{k}(\gamma(u-\Pi_{k}u))\|\,\|v_{h}\|$
and leads to

$\|\sigma_{h}^{*}-\sigma_{h}\|_{\mathbf{A}^{-1/2}}\leq\|e_{h}\|_{L}\leq\|e_{h}\|_{H}\leq
C_{8}/\beta_{0}\,\|\sigma-\sigma^{*}_{h}\|_{\mathbf{A}^{-1}}+\operatorname{\mathcal{L}ip}(\gamma)/\beta_{0}\,\|h_{\mathcal{T}}(u-\Pi_{k}u)\|.$

A triangle inequality in $L^{2}(\Omega;\mathbb{R}^{n})$ is followed by (16) in
the proof of

$\displaystyle\|\sigma-\sigma_{h}\|_{\mathbf{A}^{-1/2}}$ $\displaystyle\leq
C_{3}(1+C_{8}/\beta_{0})\Bigl{(}\operatorname{dist}(\sigma,M_{h})+\|h_{\mathcal{T}}(1-\Pi_{k})\operatorname{div}\sigma\|\Bigr{)}$
$\displaystyle\qquad+\operatorname{\mathcal{L}ip}(\gamma)/\beta_{0}\,\|h_{\mathcal{T}}(u-\Pi_{k}u)\|.$

This completes the rest of the proof. $\qquad\Box$

###  4.3 Conservative formulation

Theorem 2.a includes an error estimate for the conservative formulation with
${\mathbf{b}}_{1}:=\mathbf{A}^{-1}{\mathbf{b}}$ and ${\mathbf{b}}_{2}:=0$ in
(4) and $\sigma=-\mathbf{A}\nabla u-u\,{\mathbf{b}}$ with
$\operatorname{div}\sigma=f-\gamma u$. The refined analog of Theorem 2.b is
not expected because of an extra term exemplified in the extreme case of
piecewise constant coefficients $\mathbf{b}_{1}$ and $\gamma$. The arguments
of Subsection 4.1 lead to

$\beta_{0}\,\|e_{h}\|_{H}\leq
b((\sigma-\sigma_{h}^{*},u-\Pi_{k}u),y_{h})=(\sigma-\sigma^{*}_{h},\tau_{h})_{\mathbf{A}^{-1}}+((u-\Pi_{k}u)\mathbf{b}_{1},\tau_{h}-\Pi_{k}\tau_{h})_{L^{2}(\Omega)}.$

The last term is not of higher order for the BDM finite element family as
pointed out in [8] through numerical evidence. For the RT finite element
family, however, Lemma 2.9 shows
$\|\tau_{h}-\Pi_{k}\tau_{h}\|{\lesssim}\|h_{\mathcal{T}}\,\operatorname{div}\tau_{h}\|$
and then leads to a higher-order contribution in the asserted inequality of
Theorem 1.2.b as the final result. $\Box$

The arguments could be generalised, but those result are of limited relevance
as the convergence order is not generally improved in comparison with Theorem
2.a. An exception is the example of [7, Sect 3.5] (with $\mathbf{b}=0=\gamma$
on the unit ball) when Theorem 1.2.b guarantees $O(\delta^{2})$ for the
$L^{2}$ flux error for $k=0$.

## Acknowledgements

The research of the first author has been supported by the Deutsche
Forschungsgemeinschaft in the Priority Program 1748 under the project
”foundation and application of generalied mixed FEM towards nonlinear problems
in solid mechanics” (CA 151/22-2). The finalization of this paper has been
supported by DST SERB MATRICS grant No. MTR/2017/000199 (NN), MATRICS grant
No. MTR/201S/000309 (AKP) and SPARC project (id 235) entitled the mathematics
and computation of plates.

## References

  * [1] T. Arbogast and Z. Chen, _On the implementation of mixed methods as nonconforming methods for second-order elliptic problems_ , Math. Comp. 64 (1995), no. 211, 943–972.
  * [2] D. Boffi, F. Brezzi, and M. Fortin, _Mixed finite element methods and applications_ , Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013.
  * [3] D. Braess, _Finite elements_ , third ed., Cambridge University Press, Cambridge, 2007.
  * [4] S. C. Brenner and L. R. Scott, _The mathematical theory of finite element methods_ , third ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008.
  * [5] C. Carstensen, A. K. Dond, N. Nataraj, and A. K. Pani, _Error analysis of nonconforming and mixed FEMs for second-order linear non-selfadjoint and indefinite elliptic problems_ , Numer. Math. 133 (2016), no. 3, 557–597.
  * [6] C. Carstensen, D. Gallistl, and M. Schedensack, _$L^{2}$ best approximation of the elastic stress in the Arnold-Winther FEM_, IMA J. Numer. Anal. 36 (2016), no. 3, 1096–1119.
  * [7] C. Carstensen, D. Peterseim, and M. Schedensack, _Comparison results of finite element methods for the Poisson model problem_ , SIAM J. Numer. Anal. 50 (2012), no. 6, 2803–2823.
  * [8] A. Demlow, _Suboptimal and optimal convergence in mixed finite element methods_ , SIAM J. Numer. Anal. 39 (2002), no. 6, 1938–1953.
  * [9] J. Douglas, Jr. and J. E. Roberts, _Mixed finite element methods for second order elliptic problems_ , Mat. Apl. Comput. 1 (1982), no. 1, 91–103.
  * [10] , _Global estimates for mixed methods for second order elliptic equations_ , Math. Comp. 44 (1985), no. 169, 39–52.
  * [11] A. Ern, T. Gudi, I. Smears, and M. Vohralik, _Equivalence of local- and global-best approximations, a simple stable local commuting projector, and optimal hp approximation estimates in H(div)_ , IMA Journal of Numerical Analysis (2021).
  * [12] A. Ern and J-L . Guermond, _Theory and practices of finite elements_ , first ed., Springer-Verlag, New York, 2004.
  * [13] A. Ern and J.-L. Guermond, _Mollification in strongly Lipschitz domains with application to continuous and discrete de Rham complexes_ , Comput. Methods Appl. Math. 16 (2016), no. 1, 51–75.
  * [14] L. Gastaldi and R. H. Nochetto, _On $L^{\infty}$-accuracy of mixed finite element methods for second order elliptic problems_, Mat. Apl. Comput. 7 (1988), no. 1, 13–39.
  * [15] J. Huang and Y. Xu, _Convergence and complexity of arbitrary order adaptive mixed element methods for the Poisson equation_ , Sci. China Math. 55 (2012), no. 5, 1083–1098.
  * [16] L.E. Payne and H.F. Weinberger, _An optimal Poincaré inequality for convex domains_ , Archive for Rational Mechanics and Analysis 5 (1960), no. 1, 286–292.
  * [17] A. H. Schatz and J. P. Wang, _Some new error estimates for Ritz-Galerkin methods with minimal regularity assumptions_ , Math. Comp. 65 (1996), no. 213, 19–27.

## Appendix: Proof of Lemma 2.9

For any simplex $T\subset\mathbb{R}^{n}$, let $P_{k}(T;\mathbb{R}^{n})$ be the
linear space of vector-valued polynomials $q_{k}$ of degree at most $k$ in any
component and let $\|\bullet\|$ abbreviate the $L^{2}$ norm
$\|\bullet\|_{L^{2}(T)}$ on $T$. The particular structure of the RT function
$\tau_{RT}$ leads to some polynomial $g\in P_{k}(T)$ and

$\tau_{RT}=g(x)\,x+p_{k}\quad\text{for all }x\in T\quad\text{ and some
}p_{k}\in P_{k}(T;\mathbb{R}^{n}).$

The argument $x$ (will always belong to $T$) is often neglected as in
$\tau_{RT}:=\tau_{RT}(x)$ or $p_{k}=p_{k}(x)$, while (with a small
inconsistency, but the right emphasis) written out in the leading term
$g(x)\,x$. The latter polynomial is either identically zero or of exact degree
$k+1$ in the sense that $g$ is a sum of monomials of exact degree $k$. Adopt a
multi-index notation with
$\alpha=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}_{0}^{n}$ and the monomial
$x^{\alpha}:=x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}\cdots x_{n}^{\alpha_{n}}$ of
degree $k=|\alpha|:=\alpha_{1}+\dots+\alpha_{n}$ for any
$x=(x_{1},\dots,x_{n})\in T$. With real coefficients $c_{\alpha}$ for any
$\alpha\in\mathbb{N}_{0}^{n}$ of degree $|\alpha|=k$,

$g(x)=\sum_{|\alpha|=k}c_{\alpha}x^{\alpha}\quad\text{for all }x\in T.$

(The symbol $|\alpha|=k$ under the sum sign abbreviates the set of all multi-
indices $\alpha=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}_{0}^{n}$ of degree
$k$). The divergence

$\operatorname{div}(g(x)\,x)=n\,g(x)+x\cdot\nabla g(x)$

of the vector-valued polynomial $g(x)\,x$ of degree $k+1$ with respect to $x$
is computed with the observation that, for any $\alpha\in\mathbb{N}_{0}^{n}$
with $|\alpha|=k$,

$x\cdot\nabla(x^{\alpha})=\sum_{j=1}^{n}x_{j}\partial(x^{\alpha})/\partial
x_{j}=\sum_{j=1}^{n}\alpha_{j}x^{\alpha}=k\,x^{\alpha}.$

Consequently, $x\cdot\nabla g(x)=k\,g(x)$ and

$\operatorname{div}(g(x)\,x)=(n+k)\,g(x)\quad\text{for all }x\in T.$

This proves
$\operatorname{div}\tau_{RT}=\operatorname{div}(g(x)\,x)+q_{k-1}=(n+k)\,g(x)+q_{k-1}$
for some $q_{k-1}\in P_{k-1}(T)$. The comparison with
$\tau_{RT}=g(x)\,x+p_{k}$ leads to some polynomial remainder $r_{k}\in
P_{k}(T;\mathbb{R}^{n})$ in

$\tau_{RT}=(n+k)^{-1}\,(\operatorname{div}\tau_{RT})\,x+r_{k}\quad\text{for
all }x\in T.$

In other words, since $\operatorname{div}\tau_{RT}\in P_{k}(T)$,

$(n+k)(1-\Pi_{k})\tau_{RT}=(1-\Pi_{k})\left((\operatorname{div}\tau_{RT})\,x\right)=(1-\Pi_{k})\left((\operatorname{div}\tau_{RT})\,(x-c)\right)$

for any constant vector $c$. For instance, the center of inertia
$c=\operatorname{mid}(T)$ of $T$ with diameter $h_{T}$ satisfies
$|x-\operatorname{mid}(T)|\leq\left(n/(n+1)\right)h_{T}$ for all $x\in T$.
This leads to

$\displaystyle(n+k)\|$
$\displaystyle\tau_{RT}-\Pi_{k}\tau_{RT}\|=\|(1-\Pi_{k})\left((\operatorname{div}\tau_{RT})(x-\operatorname{mid}(T))\right)\|$
$\displaystyle\leq\|(\operatorname{div}\tau_{RT})(x-\operatorname{mid}(T))\|\leq\left(n/(n+1)\right)h_{T}\,\|\operatorname{div}\tau_{RT}\|.$

This proves $\|\tau_{RT}-\Pi_{k}\tau_{RT}\|_{L^{2}(T)}\leq\frac{n\,\
h_{T}}{(n+1)(n+k)}\|\operatorname{div}\tau_{RT}\|_{L^{2}(T)}$. $\Box$


11institutetext: 1Yale University, New Haven, CT 06511, USA
2Rensselaer Polytechnic Institute, Troy, NY 12180, USA
11email: {huidong.xie<EMAIL_ADDRESS>
# Transformer-based Dual-domain Network for Few-view Dedicated Cardiac SPECT
Image Reconstructions

Huidong Xie1 Bo Zhou1 Xiongchao Chen1 Xueqi Guo1 Stephanie Thorn1 Yi-Hwa Liu1
Ge Wang2 Albert Sinusas1 Chi Liu1

###### Abstract

Cardiovascular disease (CVD) is the leading cause of death worldwide, and
myocardial perfusion imaging using SPECT has been widely used in the diagnosis
of CVDs. The GE 530/570c dedicated cardiac SPECT scanners adopt a stationary
geometry to simultaneously acquire 19 projections to increase sensitivity and
achieve dynamic imaging. However, the limited amount of angular sampling
negatively affects image quality. Deep learning methods can be implemented to
produce higher-quality images from stationary data. This is essentially a few-
view imaging problem. In this work, we propose a novel 3D transformer-based
dual-domain network, called TIP-Net, for high-quality 3D cardiac SPECT image
reconstructions. Our method aims to first reconstruct 3D cardiac SPECT images
directly from projection data without the iterative reconstruction process by
proposing a customized projection-to-image domain transformer. Then, given its
reconstruction output and the original few-view reconstruction, we further
refine the reconstruction using an image-domain reconstruction network.
Validated by cardiac catheterization images, diagnostic interpretations from
nuclear cardiologists, and defect size quantified by an FDA 510(k)-cleared
clinical software, our method produced images with higher cardiac defect
contrast on human studies compared with previous baseline methods, potentially
enabling high-quality defect visualization using stationary few-view dedicated
cardiac SPECT scanners.

###### Keywords:

Cardiac SPECT Few-view Imaging Transformer

##  1 Introduction

The GE Discovery NM Alcyone 530c/570c [1] are dedicated cardiac SPECT systems
with 19 cadmium zinc telluride (CZT) detector modules designed for stationary
imaging. Limited amount of angular sampling on scanners of this type could
affect image quality. Due to the unique geometry of Alcyone scanners, the
centers of FOV vary at different angular positions. Hence, unlike CT scanners,
there is no straightforward method to combine projections at different
positions on Alcyone scanners. Xie et al. [12] proposed to incorporate the
displacements between centers of FOV and the rotation angles into the
iterative method for multi-angle reconstructions with registration steps
within each iteration. Despite its promising results, acquiring multi-angle
projections on this scanner is time-consuming and inconvenient in reality.
Rotating the detectors also limits the capability of dynamic imaging. Thus, it
is desirable to obtain high-quality rotation-based reconstruction directly
from the stationary SPECT projection data. This is essentially a few-view
reconstruction problem.

Previous works have attempted to address this problem by using deep-learning-
based image-to-image networks. Xie et al. proposed a 3D U-net-like network to
directly synthesize dense-view images from few-view counterparts [12]. Since
convolutional networks have limited receptive fields due to small kernel size,
Xie et al. further improved their method with a transformer-based image-to-
image network for SPECT reconstruction [13]. Despite their promising
reconstruction results, these methods use MLEM (maximum likelihood expectation
maximization) reconstruction from one-angle few-view acquisition data as
network input. The initial MLEM reconstruction may contain severe image
artifacts with important image features lost during the iterative
reconstruction process, thus would be challenging to be recovered with image-
based methods. Learning to reconstruct images directly from the projection
data could lead to improved quality.

There are a few previous studies proposed to learn the mapping between raw
data and images. AUTOMAP [14] utilized fully-connected layers to learn the
inverse Fourier transform between k-space data and the corresponding MRI
images. While such a technique could be theoretically adapted to other imaging
modalities, using a similar approach would require a significant amount of
trainable parameters, and thus is infeasible for 3D data. Würfl et al. [10]
proposed a back-projection operator to link projections and images to reduce
memory burden for CT. However, their back-projection process is not learnable.
There are also recent works [5, 6, 11] that tried to incorporate the embedded
imaging physics and geometry of CT scanners into the neural networks to reduce
redundant trainable parameters for domain mapping. However, these networks are
unable to be generalized to other imaging modalities due to different physical
properties. Moreover, these methods are hard to be extended to 3D cases
because of the geometric complexity and additional memory/computational
burden.

Here, we propose a novel 3D Transformer-based Dual-domain (projection & image)
Network (TIP-Net) to address these challenges. The proposed method
reconstructs high-quality few-view cardiac SPECT using a two-stage process.
First, we develop a 3D projection-to-image transformer reconstruction network
that directly reconstructs 3D images from the projection data. In the second
stage, this intermediate reconstruction is combined with the original few-view
reconstruction for further refinement, using an image-domain reconstruction
network. Validated on physical phantoms, porcine, and human studies acquired
on GE Alcyone 570c SPECT/CT scanners, TIP-Net demonstrated superior
performance than previous baseline methods. Validated by cardiac
catheterization (Cath) images, diagnostic results from nuclear cardiologists,
and cardiac defect quantified by an FDA-510(k)-cleared software, we also show
that TIP-Net produced images with higher resolution and cardiac defect
contrast on human studies, as compared to previous baselines. Our method could
be a clinically useful tool to improve cardiac SPECT imaging.

##  2 Methodology

Following the acquisition protocol described by Xie et al. [12], we acquired a
total of eight porcine and two physical phantom studies prospectively with
four angles of projections. The training target in this work is four-angle
data with 76 ($19\times 4$) projections and the corresponding input is one-
angle data with 19 projections. Twenty clinical anonymized 99mTc-tetrofosmin
SPECT human studies were retrospectively included for evaluation. Since only
one-angle data was available for the human studies, Cath images and clinical
interpretations from nuclear cardiologists were used to determine the
presence/absence of true cardiac defects and to assess the images
reconstructed by different methods. The use of animal and human studies was
approved by the Institutional Animal Care & Use Committee (IACUC) of Yale
University.

###  2.1 Network Structure

Figure 1: Overall TIP-Net structure. The TIP-Net is divided into P-net and
I-net. The back-projected image volume was generated by multiplying the 3D
projections with the system matrix. FC layer: fully-connected layer. Both
3D-$\mathrm{CNN_{1}}$ and 3D-$\mathrm{CNN_{2}}$ share the same structure but
with different trainable parameters.

The overall network structure is presented in Fig. 1. TIP-Net is divided into
two parts. The transformer-based [2] projection-net (P-net) aims at
reconstructing 3D SPECT volumes directly from 3D projections obtained from the
scanner. Information from the system matrix ($\mathrm{IMG_{bp}}$) is also
incorporated in P-net as prior knowledge for image reconstruction.
$\mathrm{IMG_{bp}}$ is obtained by multiplying the system matrix
$\mathrm{S}\in\mathbb{R}^{19,456\times 245,000}$ with the projection data. The
outputs from P-net ($\mathrm{IMG_{p}}$) serve as an input for image-net
(I-net).

To reconstruct $\mathrm{IMG_{p}}$ in P-net, we may simply project the output
from the Transformer block to the size of 3D reconstructed volume (i.e.,
$70\times 70\times 50$). However, such implementation requires a significant
amount of GPU memory. To alleviate the memory burden, we proposed to learn the
3D reconstruction in a slice-by-slice manner. The output from the Transformer
block is projected to a single slice ($70\times 70$), and the whole
Transformer network (red rectangular in Fig.1) is looped 50 times to produce a
whole 3D volume ($70\times 70\times 50$). Since different slice has different
detector sensitivity, all the 50 loops use different trainable parameters and
the $i^{th}$ loop aims to reconstruct the $i^{th}$ slice in the 3D volume.
Within each loop, the Transformer block takes the entire 3D projections as
input to reconstruct the $i^{th}$ slice in the SPECT volume. With the self-
attention mechanism, all 50 loops can observe the entire 3D projections for
reconstruction. The $70\times 70$ slice is then combined with
$\mathrm{IMG_{bp}}$ (only at $i^{th}$ slice), and the resized 3D projection
data. The resulting feature maps ($70\times 70\times 21$) are fed into a
shallow 2D CNN to produce a reconstructed slice at the $i^{th}$ loop
($70\times 70\times 1$), which is expected to be the $i^{th}$ slice in the
SPECT volume.

I-net takes images reconstructed by P-net ($\mathrm{IMG_{p}}$), prior
knowledge from the system matrix ($\mathrm{IMG_{bp}}$), and images
reconstructed by MLEM ($\mathrm{IMG_{mlem}}$) using one-angle data for final
reconstructions. With such a design, the network can combine information from
both domains for potentially better reconstructions compared with methods that
only consider single-domain data.

Both $\mathrm{3D\text{-}CNN_{1}}$ and $\mathrm{3D\text{-}CNN_{2}}$ use the
same network structure as the network proposed in [12], but with different
trainable parameters.

###  2.2 Optimization, Training, and Testing

In this work, the TIP-Net was trained using a Wasserstein Generative
Adversarial Network (WGAN) with gradient penalty [4]. A typical WGAN contains
2 separate networks, one generator ($G$) and a discriminator ($D$). In this
work, the network $G$ is the TIP-Net depicted in Fig.1. Formulated below, the
overall objective function for network $G$ includes SSIM, MAE, and Adversarial
loss.

$\begin{split}\underset{{\theta}_{G}}{\min}\,L=\bigg{\\{}\ell(G(I_{\mathrm{one}}),I_{\mathrm{four}})+\lambda_{a}\ell(\mathrm{P_{net}}(I_{\mathrm{one}}),I_{\mathrm{four}})-\underbrace{\lambda_{b}\,\mathbb{E}_{I_{\mathrm{one}}}\left[D(G({I}_{\mathrm{one}}))\right]}_{\mathrm{Adversarial\
Loss}}\bigg{\\}},\end{split}$ (1)

where $I_{\mathrm{one}}$ and $I_{\mathrm{four}}$ denote images reconstructed
by one-angle data and four-angle data respectively using MLEM. $\theta_{G}$
represents trainable parameters of network $G$.
$\mathrm{P_{net}}(I_{\mathrm{one}})$ represents the output of the P-net
($\mathrm{IMG_{p}}$). The function $\ell$ is formulated as below:

$\begin{split}\ell(X,Y)=\ell_{\mathrm{MAE}}(X,Y)+\lambda_{c}\,\ell_{\mathrm{SSIM}}(X,Y)+\lambda_{d}\,\ell_{\mathrm{MAE}}(\mathrm{SO}(X),\mathrm{SO}(Y)),\end{split}$
(2)

where $X$ and $Y$ represent two image volumes used for calculations.
$\ell_{\mathrm{MAE}}$ and $\ell_{\mathrm{SSIM}}$ represent MAE and SSIM loss
functions, respectively. The Sobel operator ($\mathrm{SO}$) was used to obtain
edge images, and the MAE between them was included as the loss function.
$\lambda_{a}=0.1$, $\lambda_{b}=0.005$, $\lambda_{c}=0.8$, and
$\lambda_{d}=0.1$ were fine-tuned experimentally. Network $D$ shares the same
structure as that proposed in [12].

The Adam method [8] was used for optimizations. The parameters were
initialized using the Xavier method [3]. 250 volumes of simulated 4D extended
cardiac-torso (XCAT) phantoms [9] were used for network pre-training. Leave-
one-out testing process was used to obtain testing results for all the
studies.

###  2.3 Evaluations

Reconstructed images were quantitatively evaluated using SSIM, RMSE, and PSNR.
Myocardium-to-blood pool (MBP) ratios were also included for evaluation. For
myocardial perfusion images, higher ratios are favorable and typically
represent higher image resolution. Two ablated networks were trained and used
for comparison. One network, denoted as $\mathrm{3D\text{-}CNN}$, shared the
same structure as either $\mathrm{3D\text{-}CNN_{1}}$ or
$\mathrm{3D\text{-}CNN_{2}}$ but only with $\mathrm{IMG_{mlem}}$ as input. The
other network, denoted as $\mathrm{Dual\text{-}3D\text{-}CNN}$, used the same
structure as the TIP-Net outlined in Fig.1 but without any projection-related
inputs. Compared with these two ablated networks, we could demonstrate the
effectiveness of the proposed projection-to-image module in the TIP-Net. We
also compared the TIP-Net with another transformer-based network (SSTrans-3D)
[13]. Since SSTrans-3D only considers image-domain data, comparing TIP-Net
with SSTrans-3D could demonstrate the effectiveness of the transformer-based
network for processing projection-domain information.

To compare the performance of cardiac defect quantifications, we used the FDA
510(k)-cleared Wackers-Liu Circumferential Quantification (WLCQ) software [7]
to calculate the myocardial perfusion defect size (DS). For studies without
cardiac defect, we should expect lower measured defect size as the uniformity
of myocardium improves. For studies with cardiac defects, we should expect
higher measured defect size as defect contrast improves.

Cath images and cardiac polar maps are also presented. Cath is an invasive
imaging technique used to determine the presence/absence of obstructive
lesions that results in cardiac defects. We consider Cath as the gold standard
for the defect information in human studies. The polar map is a 2D
representation of the 3D volume of the left ventricle. All the metrics were
calculated based on the entire 3D volumes. All the clinical descriptions of
cardiac defects in this paper were confirmed by nuclear cardiologists based on
SPECT images, polar maps, WLCQ quantification, and Cath images (if available).

##  3 Results

###  3.1 Porcine results

Results from one sample porcine study are presented in Fig. S1. This pig had a
large post-occlusion defect created by inflating an angioplasty balloon in the
left anterior descending artery. As pointed out by the green arrows in Fig.
S1, all neural networks improved the defect contrast with higher MBP ratios
compared with the one-angle images. The TIP-Net results were better than other
network results in terms of defect contrast based on defect size measured by
WLCQ. TIP-Net also produced images with clearer right ventricles, as
demonstrated by the Full-width at Half Maximum (FHMW) values presented in Fig.
S1.

Quantitative results for 8 porcine and 2 physical phantom studies are included
in Table 1. Based on paired t-tests, all the network results are statistically
better than one-angle results ($p<0.001$), and the TIP-Net results are
statistically better than all the other three network results ($p<0.05$). The
average MBP ratios shown in Table 1 indicate that the proposed TIP-Net
produced images with higher resolution compared with image-based networks.

Average defect size measurements for porcine and physical phantom studies with
cardiac defects are 35.9%, 42.5%, 42.3%, 43.6%, 47.0%, 46.2% for one-angle,
3D-CNN, Dual-3D CNN, TIP-Net, and four-angle image volumes, respectively.
Images reconstructed by TIP-Net present overall higher defect contrast and the
measured defect size is closest to the four-angle images.

For normal porcine and physical phantom studies without cardiac defect, these
numbers are 16.0%, 12.0%, 14.7%, 11.2%, 11.5%, and 10.1%. All neural networks
showed improved uniformity in the myocardium with lower measured defect size.
Both transformer-based methods, TIP-Net and SSTrans-3D, showed better defect
quantification than other methods on these normal subjects.

Table 1: Quantitative values for porcine and phantom results obtained using
different methods ($\mathrm{MEAN}\pm\mathrm{STD}$). Images reconstructed using
four-angle data were used as the reference. Best values are marked in bold.
$p<0.05$ observed between all groups.

Evaluation | One-angle | 3D-CNN | Dual-3D-CNN | SSTrans-3D | TIP-Net | Four-angle
---|---|---|---|---|---|---
SSIM$\uparrow$ | $0.945\pm 0.009$ | $0.951\pm 0.008$ | $0.950\pm 0.008$ | $0.947\pm 0.009$ | $0.953\pm 0.008$ |
PSNR$\uparrow$ | $35.764\pm 3.546$ | $36.763\pm 3.300$ | $36.787\pm 3.208$ | $36.613\pm 3.509$ | $38.114\pm 2.876$ |
RMSE$\downarrow$ | $0.020\pm 0.005$ | $0.0181\pm 0.005$ | $0.0185\pm 0.005$ | $0.0194\pm 0.005$ | $0.0179\pm 0.005$ |
MBP$\uparrow$ | $5.16\pm 1.83$ | $10.73\pm 2.98$ | $10.21\pm 3.17$ | $11.10\pm 2.78$ | $12.73\pm 3.73$ |  $13.63\pm 3.91$

###  3.2 Human results

Chosen by a nuclear cardiologist, results from one representative human study
with multiple cardiac defects are presented in Fig. 2. Note that there was no
four-angle data for human studies. Based on clinical interpretations from the
nuclear cardiologist and Cath results, this patient had multiple perfusion
defects in the apical (blue arrows in Fig. 2) and basal (green and yellow
arrows in Fig. 2) regions. As presented in Fig. 2, all the networks produced
images with better defect contrasts. The proposed TIP-Net further improved the
defect contrast, which is favorable and could make the clinical
interpretations more confident.

Figure 2: Results from a sample stress human study. Red arrows point to
different coronary vessels. Non-red arrows with the same color point to the
same defect/stenosis in the heart. Numbers in the parentheses are the MBP
ratios and the defect size (%LV).

Another human study was selected and presented in Fig. 3. This patient had
stenosis in the distal left anterior descending coronary artery, leading to a
medium-sized defect in the entire apical region (light-green arrows). As
validated by the Cath images, polar map, and interpretation from a nuclear
cardiologist, TIP-Net produced images with higher apical defect contrast
compared with other methods, potentially leading to a more confident
diagnostic decision.

Figure 3: Results from a sample stress human study. Red arrows point to
different coronary vessels. Light-green arrows point to the same
defect/stenosis in the heart. Numbers in the parentheses are the MBP ratios
and the defect size (%LV).

The average MBP ratios on human studies for one-angle, 3D-CNN, Dual-3D-CNN,
SSTrans-3D, and TIP-Net images are $3.13\pm 0.62$, $4.08\pm 0.83$, $4.31\pm
1.07$, $4.17\pm 0.99$, and $4.62\pm 1.29$, respectively
($\mathrm{MEAN}\pm\mathrm{STD}$). The higher ratios in images produced by TIP-
Net typically indicate higher image resolution.

14 patients in the testing data have cardiac defects, according to diagnostic
results. The average defect size measurements for these patients are 16.8%,
22.6%, 21.0%, 22.4%, and 23.6% for one-angle, 3D-CNN, Dual-3D-CNN, SSTrans-3D,
and TIP-Net results. The higher measured defect size of the TIP-Net indicates
that the proposed TIP-Net produced images with higher defect contrast.

For the other 6 patients without cardiac defects, these numbers are 11.5%,
12.8%, 13.8%, 12.4%, and 11.8%. These numbers show that TIP-Net did not
introduce undesirable noise in the myocardium, maintaining the overall
uniformity for normal patients. However, other deep learning methods tended to
introduce non-uniformity in these normal patients and increased the defect
size.

###  3.3 Intermediate network output

To further show the effectiveness of P-net in the overall TIP-Net design,
outputs from P-net (images reconstructed by the network directly from
projections) are presented in Fig. S2. The presented results demonstrate that
the network can learn the mapping between 3D projections to 3D image volumes
directly without the iterative reconstruction process. In the porcine study,
limited angular sampling introduced streak artifacts in the MLEM-reconstructed
one-angle images (yellow arrows in Fig. S2). P-net produced images with fewer
few-view artifacts and higher image resolution.

The human study presented in Fig. S2 has an apical defect according to the
diagnostic results (blue arrows in Fig. S2). However, this apical defect is
barely visible in the one-angle image. P-net produced images with higher
resolution and improved defect contrast. Combining both outputs (one-angle
images and $\mathrm{IMG_{p}}$), TIP-Net further enhanced the defect contrast.

##  4 Discussion and Conclusion

We proposed a novel TIP-Net for 3D cardiac SPECT reconstruction. To the best
of our knowledge, this work is the first attempt to learn the mapping from 3D
realistic projections to 3D image volumes. Previous works in this direction
[14, 6, 5] were 2D and only considered simulated projections with ideal
conditions. 3D realistic projection data have more complex geometry are also
affected by other physical factors that may not exist in simulated
projections.

The proposed method was tested for myocardial perfusion SPECT imaging.
Validated by nuclear cardiologists, diagnostic results, Cath images, and
defect size measured by WLCQ, the proposed TIP-Net produced images with higher
resolution and higher defect contrast for patients with perfusion defects. For
normal patients without perfusion defects, TIP-Net maintained overall
uniformity in the myocardium with higher image resolution. Similar performance
was observed in porcine and physical phantom studies.

## References

  * [1] Bocher, M., Blevis, I., Tsukerman, L., Shrem, Y., Kovalski, G., Volokh, L.: A fast cardiac gamma camera with dynamic SPECT capabilities: design, system validation and future potential. Eur J Nucl Med Mol Imaging 37(10), 1887–1902 (Oct 2010), https://doi.org/10.1007/s00259-010-1488-z
  * [2] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., Houlsby, N.: An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale. arXiv:2010.11929 [cs] (Jun 2021), http://arxiv.org/abs/2010.11929
  * [3] Glorot, X., Bengio, Y.: Understanding the difficulty of training deep feedforward neural networks. In: Teh, Y.W., Titterington, M. (eds.) Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research, vol. 9, pp. 249–256. PMLR, Chia Laguna Resort, Sardinia, Italy (13–15 May 2010), http://proceedings.mlr.press/v9/glorot10a.html
  * [4] Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., Courville, A.: Improved training of wasserstein gans. In: Advances in Neural Information Processing Systems. vol. 30. Curran Associates, Inc. (2017)
  * [5] He, J., Wang, Y., Ma, J.: Radon Inversion via Deep Learning. IEEE Transactions on Medical Imaging 39(6), 2076–2087 (Jun 2020). https://doi.org/10.1109/TMI.2020.2964266
  * [6] Li, Y., Li, K., Zhang, C., Montoya, J., Chen, G.: Learning to Reconstruct Computed Tomography (CT) Images Directly from Sinogram Data under A Variety of Data Acquisition Conditions. IEEE Transactions on Medical Imaging pp. 1–1 (2019). https://doi.org/10.1109/TMI.2019.2910760
  * [7] Liu, Y., Sinusas, A., DeMan, P., Zaret, B., Wackers, F.: Quantification of SPECT myocardial perfusion images: methodology and validation of the Yale-CQ method. J Nucl Cardiol 6(2), 190–204 (Apr 1999). https://doi.org/10.1016/s1071-3581(99)90080-6
  * [8] PK, D., B, J.: Adam: A method for stochastic optimization. In: 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (2015), http://arxiv.org/abs/1412.6980
  * [9] Segars, W., Sturgeon, G., Mendonca, S., Grimes, J., Tsui, B.: 4D XCAT phantom for multimodality imaging research. Medical Physics 37(9), 4902–4915 (2010). https://doi.org/10.1118/1.3480985, https://aapm.onlinelibrary.wiley.com/doi/abs/10.1118/1.3480985
  * [10] Würfl, T., Hoffmann, M., Christlein, V., Breininger, K., Huang, Y., Unberath, M., Maier, A.K.: Deep Learning Computed Tomography: Learning Projection-Domain Weights From Image Domain in Limited Angle Problems. IEEE Transactions on Medical Imaging 37(6), 1454–1463 (Jun 2018). https://doi.org/10.1109/TMI.2018.2833499
  * [11] Xie, H., Shan, H., Cong, W., Liu, C., Zhang, X., Liu, S., Ning, R., Wang, G.: Deep Efficient End-to-End Reconstruction (DEER) Network for Few-View Breast CT Image Reconstruction. IEEE Access 8, 196633–196646 (2020). https://doi.org/10.1109/ACCESS.2020.3033795
  * [12] Xie, H., Thorn, S., Chen, X., Zhou, B., Liu, H., Liu, Z., Lee, S., Wang, G., Liu, Y.H., Sinusas, A.J., Liu, C.: Increasing angular sampling through deep learning for stationary cardiac SPECT image reconstruction. Journal of Nuclear Cardiology 30(1), 86–100 (2022). https://doi.org/10.1007/s12350-022-02972-z
  * [13] Xie, H., Thorn, S., Liu, Y.H., Lee, S., Liu, Z., Wang, G., Sinusas, A.J., Liu, C.: Deep learning based few-angle cardiac SPECT reconstruction using transformer. IEEE Transactions on Radiation and Plasma Medical Sciences 7(1), 33–40 (2022). https://doi.org/10.1109/TRPMS.2022.3187595
  * [14] Zhu, B., Liu, J.Z., Cauley, S.F., Rosen, B.R., Rosen, M.S.: Image reconstruction by domain-transform manifold learning. Nature 555(7697), 487–492 (2018). https://doi.org/10.1038/nature25988

Supplemental Materials

##  5 Additional Figures

Figure S1: A sample porcine study. Green arrows point to the post-occlusion
reperfusion defect in the heart. Numbers in the parentheses are the MBP ratios
and the defect size (%LV). FWHM values were calculated along the dashed blue
lines through the RV myocardium in SA slices. VLA: vertical-long axis. SA:
short-axis.

Figure S2: Network outputs at different stages in TIP-Net. Yellow arrows in
pig study point to streak artifacts were introduced due to limited angular
sampling in one-angle images. Blue arrows point to the cardiac defect,
confirmed by cardiologist, which is barely invisible in the one-angle image.
Numbers in the parentheses are the calculated MBP ratios and the defect size
(%LV). Defect size measurements were not available for intermediate network
output (i.e., $\mathrm{IMG_{p}}$).