\documentclass[reqno]{amsart}
\usepackage{tikz}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
Tenth MSU Conference on Differential Equations and Computational
Simulations. \newline
\emph{Electronic Journal of Differential Equations},
Conference 23 (2016),  pp. 59--76.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University.}
\vspace{9mm}}

\begin{document} \setcounter{page}{59}
\title[\hfilneg EJDE-2016/Conf/23 \hfil 
 Distributions and the Stability of DG Methods]
{Distributional derivatives and  stability of discontinuous 
 Galerkin finite element approximation methods}


\author[T. Lewis \hfil EJDE-2016/conf/23 \hfilneg]
{Thomas Lewis} 

\dedicatory{Dedicated to Professor Ratnasingham Shivaji on  his 60-th birthday}


\address{Thomas Lewis \newline
Department of Mathematics and Statistics,
The University of North Carolina at Greensboro,
 Greensboro, NC 27412, USA}
\email{tllewis3@uncg.edu}

\thanks{Published March 21, 2016.}
\subjclass[2010]{65J02, 65N02, 35C02}
\keywords{Discontinuous Galerkin methods; broken Sobolev spaces;
\hfill\break\indent  distributional derivatives; numerical derivatives}

\begin{abstract}
 The goal of this article is to explore and motivate stabilization requirements
 for various types of discontinuous Galerkin (DG) methods.
 A new approach for the understanding of DG  approximation methods for second 
 order elliptic partial differential equations is introduced.
 The approach explains the weaker stability requirements for local 
 discontinuous Galerkin (LDG) methods  when compared to interior-penalty 
 discontinuous Galerkin methods  while also motivating the existence of methods
 such as the minimal dissipation LDG method  that are stable without the addition
 of interior penalization.  The main idea is to relate the underlying DG gradient 
 approximation to   distributional  derivatives instead of the traditional 
 piecewise gradient operator associated with  broken Sobolev spaces.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks

\newcommand{\tbar}[1]{|\hspace{-0.045cm}|\hspace{-0.045cm}|{#1}|
\hspace{-0.045cm}|\hspace{-0.045cm}|}


\section{Introduction} \label{intro_sec}

In this article, a new approach for understanding discontinuous Galerkin (DG) 
finite element approximation methods for elliptic partial differential equations 
(PDEs) will be motivated and explored with regards to stability requirements.
The main emphasis will be on the difficulties associated with extending weak 
solution theory for PDEs to broken Sobolev spaces.
To this end, existence and uniqueness results, as well as Friedrichs inequalities,
 will be established for second order elliptic problems posed in broken Sobolev spaces.
Once the extension has been made, we will be able to directly recover many existing
discontinuous Galerkin approximation methods by using simple projection operators.
We will gain new insight into the increased necessity for stabilization using 
penalization techniques for Interior-Penalty Discontinuous Galerkin (IPDG) methods,
\cite{DouglasDupont76,Wheeler78,Baker77,RiviereWheelerGirault00,RiviereWheelerGirault01},
versus the weaker penalization requirements for
Local Discontinuous Galerkin (LDG) methods, \cite{Cockburn98,Castilo00}.
The approach will also motivate the existence of penalty-free methods such as
the Minimal Dissipation Local Discontinuous Galerkin (MD-LDG) method, 
\cite{CockburnDong07},
and the Dual-Wind Discontinuous Galerkin (DWDG) method,
 \cite{LewisNeilan12,feng_lewis_wise15}.
This paper is meant as a survey paper that can shed new light on known 
properties of DG methods for researchers familiar with the topic
while also motivating DG approximation methods for researchers familiar
 with elliptic PDE theory.

Let $H^m(\Omega)$ denote the set of all functions in $L^2(\Omega)$ whose 
distributional derivatives up to order $m$ are in $L^2(\Omega)$.
Define $L^2_g(\Omega)$ as the set of all functions in $L^2(\Omega)$ 
that have (a well-defined) trace value $g$ for $g \in L^2(\partial \Omega)$
and let $H^m_0(\Omega)$ denote the set of all functions in $H^m(\Omega)$ whose 
traces vanish up to order ($m-1$) on $\partial \Omega$.
For transparency, we will refer to Poisson's equation with Dirichlet boundary 
conditions as our prototypical elliptic boundary value problem:
\begin{subequations} \label{BVP}
\begin{gather}
	- \Delta u  = f ,  \quad \text{in } \Omega,  \\
	u = g ,  \quad \text{on } \partial \Omega,
\end{gather}
\end{subequations}
where $\Omega \subset \mathbb{R}^d$, $d \geq 2$, is an open, convex, 
polygonal domain;
$g \in H^{1/2}(\partial \Omega)$;
and $\Delta \equiv \sum_{j=1}^d \frac{\partial^2}{\partial {x_j}^2}$
denotes the Laplacian operator.
In weak form, \eqref{BVP} is equivalent to finding the unique function
$u \in H^1(\Omega) \cap L^2_g(\partial \Omega)$ such that
\begin{equation} \label{BVP_weak}
	\big( \nabla u , \nabla \varphi \big)_{\Omega}
= \big( f , \varphi \big)_{\Omega}
	\quad \forall \varphi \in H^1_0(\Omega) ,
\end{equation}
where $( \cdot , \cdot )_{S}$ denotes the $L^2$ inner product on $S$.

To approximate the solution to \eqref{BVP_weak} using discontinuous Galerkin
methods, we first need to formalize the problem using broken Sobolev spaces.
Traditionally, this can be done using a piecewise gradient operator
and controlling the jumps/discontinuities along interior faces/edges.
Such a framework naturally leads to IPDG methods,
and, using the unified framework in \cite{unified},
can be extended to LDG methods.
In this paper, we consider a different approach that is based on using
distributional derivatives instead of the piecewise gradient operator.
Following the example of the unified framework,
we will consider the flux-based formulation of LDG-based methods
and the Bassi-Rebay method, \cite{bassi_rebay},
where the derivative approximation is based on a well-chosen single-valued trace.
However, we will formally consider the consequences of such a flux-based
definition to build intuition and then rigorously consider the relationship of such a flux-based
approach within the setting of distributional derivatives.
The approach in this paper
complements the Hybrid Discontinuous Galerkin method, \cite{cockburn2009unified},
where the interior fluxes are directly approximated and used to determine the approximate gradient.
We will see that the new interpretation of differential operators for broken Sobolev functions
inherently motivates how to control
the jumps/discontinuities along interior faces/edges.

The remainder of this article is organized as follows.
In Section~\ref{broken_spaces} we will introduce broken Sobolev spaces
and some known results.
The IPDG method will be introduced in Section~\ref{IPDG_sec}
as well as several applications to broken Sobolev spaces.
The section will provide an introduction to DG methods
and further develop the background for broken Sobolev spaces.
A heuristic framework will be considered in Section~\ref{heuristic_sec}
that will motivate the analysis in Section~\ref{LDG}
where we consider distributional derivatives and how they relate to the LDG method.
We will derive a stability result for distributional derivatives with regard 
to broken Sobolev spaces
and then use projection operators to formalize DG methods.
We end the paper with some concluding remarks in Section~\ref{conclusion_sec}.
Note that the ideas presented in this paper are readily extended
to a general class of second order elliptic PDEs,
and the ideas
formally motivate both the Weak Galerkin Finite Element Method introduced 
in \cite{wang2013weak}
and DG methods that can be expressed using the DG finite element 
differential calculus proposed in \cite{FengLewisNeilan13}.

\section{Broken Sobolev spaces} \label{broken_spaces}

We now recall the traditional extension of \eqref{BVP_weak} to broken Sobolev spaces.
First, we must introduce some (standard) definitions and notation.
Then we will pose \eqref{BVP} in the setting of broken Sobolev spaces
and discuss the major obstacle that must be overcome.
The fundamental result will be an extension of Friedrichs inequality 
to broken Sobolev spaces, a result that will be reinterpreted in the 
context of distributional derivatives in Section~\ref{LDG}
where we successfully extend \eqref{BVP_weak} to broken Sobolev spaces
while preserving existence and uniqueness properties.

\subsection{Definitions and notation}

Let $\mathcal{T}_h$ denote a locally quasi-uniform and shape-regular simplicial triangulation of $\Omega$
(see \cite{Ciarlet78,Brenner}).
Let $\mathcal{E}_h$ denote the set of all ($d$-1)-dimensional simplices in the triangulation,
$\mathcal{E}_h^I \subset \mathcal{E}_h$ denote the set of all interior ($d$-1)-dimensional simplices in the triangulation, and
$\mathcal{E}_h^B \subset \mathcal{E}_h$ denote the set of all boundary ($d$-1)-dimensional simplices in the triangulation.
Notationally, we set $h_K$ as the diameter of the simplex $K \in \mathcal{T}_h$, $h_e$ the diameter of the
simplex $e \in \mathcal{E}_h$, and $h = \max_{K \in \mathcal{T}_h} h_K$.

We define the piecewise $L^2$ inner product with respect to the triangulation by
\[
	\big( v , w \big)_{\mathcal{T}_h} \equiv \sum_{K \in \mathcal{T}_h} \int_K v \, w \, dx ,
\]
and the piecewise $L^2$ inner product with respect to the collection of 
faces/edges $\mathcal{E}_h$ by
\[
	\langle v , w \rangle_{\mathcal{S}_h} \equiv \sum_{e \in \mathcal{S}_h} \int_e v \, w \, ds
\]
for all $\mathcal{S}_h \subseteq \mathcal{E}_h$.
We use bold-faced formatting to denote a vector-valued space;
for example, we have $\mathbf{L}^2(\mathcal{T}_h) \equiv \big[ L^2(\mathcal{T}_h) \big]^d$.

We define our broken Sobolev spaces with respect to the triangulation as
\[
	W^{m,p}(\mathcal{T}_h) \equiv \prod_{K \in \mathcal{T}_h} W^{m,p}(K)
\]
for all integers $m \geq 0$ and real numbers $p \in [1,\infty]$,
where
\[
	W^{m,p}(K) \equiv \left\{ v \in L^{p}(K) \mid
	D^\alpha v \in L^p(K) \text{ for all multi-indices } \alpha \text{ with } | \alpha | \leq m \right\} .
\]
The broken Sobolev space $W^{m,p}(\mathcal{T}_h)$ is a Banach space with respect to the norm
\[
	\| v \|_{W^{m,p}(\mathcal{T}_h)} \equiv
	\sum_{\substack{k=1 \\ | \alpha | \leq k}}^m 
\| D^\alpha v \|_{L^p(\mathcal{T}_h)} .
\]
We denote the Hilbert space $W^{m,2}(\mathcal{T}_h)$ by $H^m(\mathcal{T}_h)$, where the inner product
 is defined by
\[
	( v , w )_{H^m(\mathcal{T}_h)} \equiv \sum_{\substack{k=1 \\ | \alpha | \leq k}}^m
		\big( D^\alpha v , D^\alpha w \big)_{\mathcal{T}_h} .
\]
The space $W^{m,p}_0(\mathcal{T}_h)$ denotes the set of all functions in $W^{m,p}(\mathcal{T}_h)$
whose trace on $\partial \Omega$ vanishes up to order ($m$-1).
We also let $C^{m}_c(\Omega) \subset W^{m,p}_0(\mathcal{T}_h)$
denote the set of all functions in $C^m(\Omega)$ with compact support,
where $C^m(\Omega)$ denotes the set of all continuously differentiable functions
up to order $m$.

Broken Sobolev functions inherently are allowed to be multi-valued along
interior faces/edges,
provided trace values exists on each simplex in the triangulation.
Let $K^+, K^- \in \mathcal{T}_h$ and $e = \partial K^- \cap \partial K^+$.
Without a loss of generality, we assume the global labeling number of $K^+$ is smaller than
that of $K^-$.
We then define the sided-flux values for $v$ as
\begin{align*}
v^{+} \big|_e \equiv v \big|_{e \cap \partial K^+},  \quad
v^{-} \big|_e \equiv v \big|_{e \cap \partial K^-},
\end{align*}
where $v \big|_{\partial K}$ is understood to be the trace of $v$ defined on $\overline{K}$.
Suppose $K$ is a boundary simplex.
We extend the sided-flux definitions to the boundary of $\Omega$ by
\[
v^\pm \big|_{\partial K \cap \partial \Omega} \equiv v \big|_{\partial K \cap \partial \Omega} .
\]
The standard jump and average operators on $\mathcal{E}_h$ are defined by
\[
	[v]  \equiv
	\begin{cases}
	v^- - v^+ & \text{on } \mathcal{E}_h^I \\
	v & \text{on } \mathcal{E}_h^B,
	\end{cases}
	 \quad
	\{v\} \equiv \frac{v^- + v^+}{2} ,
\]
respectively.
We impose the convention that the outward normal vector on $e$,
denoted by $\mathbf{n}$,
is always given by the outward normal vector for $K^-$.

\begin{remark} \rm
We have $W^{m,p}(\Omega) \subset W^{m,p}(\mathcal{T}_h)$, where functions in $W^{m,p}(\mathcal{T}_h)$
are allowed to have discontinuities along interior faces/edges.
Thus, the underlying gradient operator for broken Sobolev functions, denoted
 $\nabla_h$, is understood piecewise with respect to the triangulation so that
 when acting on $W^{m,p}(\Omega)$ we have $\nabla_h = \nabla$.
\end{remark}

\subsection{Friedrichs inequality and interior penalization}

Using the above notation, we are ready to derive a ``weak formulation" 
for \eqref{BVP} in the broken Sobolev space $H^2(\mathcal{T}_h)$.
Let $v \in H^2(\mathcal{T}_h) \cap L^2_g(\Omega)$ such that
$-\big( \Delta v , \varphi \big)_{\mathcal{T}_h} = \big(f , \varphi \big)_{\mathcal{T}_h}$ 
for all $\varphi \in H^1_0(\mathcal{T}_h)$.
Then, we have
\begin{equation}  \label{BVP_weak_grad}
\begin{aligned}
	\big(f , \varphi \big)_{\mathcal{T}_h}
	& = - \big( \Delta v , \varphi \big)_{\mathcal{T}_h} \\
	& = \big( \nabla v , \nabla \varphi \big)_{\mathcal{T}_h}
		- \big\langle [\nabla v] , \big\{ \varphi \big\} \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
		- \big\langle \big\{ \nabla v \big\} , \big[ \varphi \big] \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
\end{aligned}
\end{equation}
for all $\varphi \in H^1_0(\mathcal{T}_h)$.
Suppose $f = 0$ and $g = 0$.
Then, $u = 0$ is the unique solution to \eqref{BVP_weak}.
In contrast, there exists infinitely many piecewise constant functions
$v \in H^2(\mathcal{T}_h) \cap L^2_0(\Omega)$ that solve
\eqref{BVP_weak_grad}.
One such example is the function
\[
	v(x) = \begin{cases}
	0 & \text{if } x \text{ is in a boundary simplex} ,  \\
	1 & \text{if } x \text{ is in an interior simplex} .
	\end{cases}
\]

\begin{remark} \rm \quad

(a)
The formulation \eqref{BVP_weak_grad}
is equivalent to the problem: find $v \in H^2(\mathcal{T}_h) \cap L^2_g(\Omega)$ such that
\begin{equation} \label{BVP_weak_grad_local}
	\big( f , \varphi \big)_K
	= - \big( \Delta v , \varphi \big)_K
	= \big( \nabla v , \nabla \varphi \big)_K
		- \big\langle \nabla v , \varphi \, \vec{n}_K \big\rangle_{\partial K}
\end{equation}
for all $\varphi \in H^1(K)$ for all $K \in \mathcal{T}_h$,
where $\vec{n}_K$ denotes the unit outward normal vector for $\overline{K}$.
Problem \eqref{BVP_weak_grad} is obtained by summing \eqref{BVP_weak_grad_local}
over all simplices $K$.

(b)
The solution is assumed to be in the space $H^2(\mathcal{T}_h)$ instead of $H^1(\mathcal{T}_h)$
in the above weak formulation.
The interior trace values of $\nabla v \cdot \mathbf{n}$ are not defined for all functions in $H^1(\mathcal{T}_h)$,
see \cite{grisvard}.
\end{remark}

The main issue in the previous example
is that \eqref{BVP_weak_grad} has no mechanism to address the potential
for discontinuities that characterizes functions in $H^2(\mathcal{T}_h)$.
Neighboring simplices in $\mathcal{T}_h$ have no communication across their shared boundary
 interface.
As a consequence, when the gradient operator is utilized in an entirely local, 
piecewise fashion with respect to $\mathcal{T}_h$,
Poisson's problem does not have a unique solution in the broken Sobolev space
$H^2(\mathcal{T}_h) \cap L_g^2(\Omega)$ for any given combination of boundary data $g$
and source data $f$.
This observation can be explained by the extension of Friedrichs inequality
to broken Sobolev spaces:

\begin{lemma}[{\cite[Lemma 2.1]{unified}}] \label{Friedrichs_standard}
There exists a positive constant $C$ independent of $h$ such that
\begin{equation} \label{poinc_standard}
	\| v \|_{L^2(\Omega)}^2 \leq C \Big(  \big\| \nabla v \big\|_{L^2(\mathcal{T}_h)}^2
		+ \sum_{e \in \mathcal{E}_h^B} h_e^{-1} \| v \|_{L^2(e)}^2
		+ \sum_{e \in \mathcal{E}_h^I} h_e^{-1} \| [v] \|_{L^2(e)}^2
		\Big)
\end{equation}
for all $v \in H^1(\mathcal{T}_h)$.
\end{lemma}

Thus, the $L^2$ norm of $\nabla_h v$ is not sufficient to control the $L^2$ norm 
of $v$ for all functions $v \in H^1(\mathcal{T}_h)$.
In order to control the $L^2$ norm of a broken Sobolev function using the $L^2$ norm of the
piecewise gradient,
we also need to control the $L^2$ norm of the interior jumps of the function.
Such an idea is the basis of IPDG methods as will be seen in the following section.
A more detailed account of Poincar\`e and Friedrichs inequalities for broken Sobolev 
spaces can be found in \cite{brenner03}.


\section{Interior-penalty discontinuous Galerkin methods} \label{IPDG_sec}

In this section, we introduce the (symmetric) IPDG method
and consider its application to broken Sobolev spaces.
The main idea for DG methods is to approximate the solution to \eqref{BVP}
with a function defined in a finite dimensional subspace of $H^2(\mathcal{T}_h)$.

\subsection{Notation}

For a fixed integer $r \geq 0$, define the standard DG
finite element space $V_{h,r} \subset C^\infty(\mathcal{T}_h) \subset H^1(\mathcal{T}_h)$ by
\[
	V_{h,r} \equiv \prod_{K \in \mathcal{T}_h} \mathbb{P}_r(K) ,
\]
where $\mathbb{P}_r(K)$ denotes the set of all polynomials on $K$ with degree 
not exceeding $r$.
The analogous vector-valued DG space is given by $\mathbf{V}_{h,r} \equiv \big[ V_{h,r} \big]^d$.
Most of the approximation methods below assume $r \geq 1$, although many DG
methods have been extended to the case $r=0$.

The key component for extending continuous results that hold in broken Sobolev spaces
to discrete spaces for DG methods
is a projection operator from $L^2(\Omega)$ to $V_{h,r}$.
A natural choice is $L^2$ projection.
We denote by $\mathcal{P}_{h,r} : L^2(\Omega) \to V_{h,r}$
the $L^2$ projection operator onto $V_{h,r}$,
which is defined by
\[
	\big( \mathcal{P}_{h,r} v , \varphi_h \big)_{\mathcal{T}_h} = \big( v , \varphi_h \big)_{\mathcal{T}_h}
	\quad \forall \varphi_h \in V_{h,r}
\]
for all $v \in L^2(\Omega)$.
We also let $\mathcal{\vec{P}}_{h,r} : \mathbf{L}^2(\Omega) \to \mathbf{V}_{h,r}$
denote the $L^2$ projection operator onto $\mathbf{V}_{h,r}$.

\subsection{Application to broken Sobolev spaces}

We now attempt to extend \eqref{BVP} to the broken Sobolev space $H^2(\mathcal{T}_h)$
using ideas from the (symmetric) IPDG method.
While the extension may still lack uniqueness when posed in the full broken 
Sobolev space $H^2(\mathcal{T}_h)$,
the extension does guarantee uniqueness
when restricting the solution and test functions to (finite dimensional)
DG finite element spaces.
Thus, the formulation is rich enough to construct sequences of functions
in broken Sobolev spaces such that the sequences converge
to the solution of the original problem \eqref{BVP}.
The major difference in the interior-penalty based formulation
and the first attempt in Section~\ref{broken_spaces}
is the introduction of interior jump stabilization terms
in the ``weak representation" of \eqref{BVP} which
allows for the application of Lemma~\ref{Friedrichs_standard}.

Define the (symmetric) bilinear form $B_h : H^2(\mathcal{T}_h) \times H^2(\mathcal{T}_h) \to \mathbb{R}$ 
by
\begin{align*}
	B_h(v,w)
	& \equiv \big(\nabla v , \nabla w \big)_{\mathcal{T}_h}
	- \big\langle \big\{ \nabla v \big\} , [w] \mathbf{n} \big\rangle_{\mathcal{E}_h}
	- \big\langle [v] , \{\nabla w\}
\cdot \mathbf{n} \big\rangle_{\mathcal{E}_h} \\
	& \quad
	+ \gamma \sum_{e \in \mathcal{E}_h} h_e^{-1} \big\langle [v] ,
 [w] \big\rangle_{e}
\end{align*}
for all $v,w \in H^2(\mathcal{T}_h)$
for some constant $\gamma > 0$ called the penalty parameter.
Define an auxiliary norm $\tbar{\cdot} : H^2(\mathcal{T}_h) \to [0,\infty)$ by
\[
	\tbar{v}^2 \equiv \sum_{K \in \mathcal{T}_h} \| \nabla v \|_K^2
	+ \sum_{e \in \mathcal{E}_h} \Big( 2 \frac{\gamma}{h_e} \| [v]
\|_{L^2(e)}^2
		+ \frac{h_e}{\gamma} \| \big\{ \nabla v \cdot \mathbf{n} \big\}
\|_{L^2(e)}^2 \Big)
\]
for all $v \in H^2(\mathcal{T}_h)$,
and, lastly,
define the linear functional $F_h : H^2(\mathcal{T}_h) \to \mathbb{R}$ by
\[
	F_h(w) \equiv
	\big( f , w \big)_{\mathcal{T}_h}
	- \big\langle g , \nabla w \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B} +
	\gamma \sum_{e \in \mathcal{E}_h^B} h_e^{-1} \big\langle g , w \big\rangle_{e}
\]
for all $w \in H^2(\mathcal{T}_h)$.
Notice that the above traces involving both function values and gradient values
are well-defined for all functions in $H^2(\mathcal{T}_h)$.
The (symmetric) IPDG method is defined by finding $u_h \in V_{h,r}$ that solves 
the finite dimensional problem
$B_h(u_h , \varphi_h) = F_h (\varphi_h)$ for all $\varphi_h \in V_{h,r}$.
Then, we have $u_h$ approximates the solution to \eqref{BVP_weak}.

We can extend \eqref{BVP} to broken Sobolev spaces by considering the 
(continuous) problem:
Find $u \in H^2(\mathcal{T}_h)$ such that $B_h (u , \varphi) = F_h (\varphi)$
for all $\varphi \in H^2(\mathcal{T}_h)$.
Notice that the Dirichlet boundary condition is enforced using a penalization 
along the boundary.
We could also consider the problem of finding $u \in H^2(\mathcal{T}_h) \cap L^2_g(\Omega)$
with the test functions $\varphi \in H^2(\mathcal{T}_h) \cap L^2_0(\Omega)$.
In this case, some of the boundary face/edge terms would disappear in the 
definitions of $B_h$, $\tbar{\cdot}$, and $F_h$.
However, when implementing the IPDG method using piecewise polynomials,
the boundary conditions are typically enforced weakly using boundary penalization
due to the fact the Dirichlet data cannot be satisfied when $g$ is not a polynomial
with degree less than or equal to $r$.

The following facts are well-documented in the DG literature:

\begin{lemma} \quad
\begin{itemize}
\item[(i)]
$\tbar{\cdot}$ is a norm on $H^2(\mathcal{T}_h)$.
\item[(ii)]
Consistency:
If $u \in H^2(\Omega) \cap L^2_g(\Omega)$ is a solution to \eqref{BVP},
then $B_h(u, \varphi) = F_h(\varphi)$ for all $\varphi \in H^2(\mathcal{T}_h)$.

\item[(iii)]
Continuity:
There exists a positive constant $C_{B,1}$ such that, for all $v,w \in H^2(\mathcal{T}_h)$,
there holds $\left| B_h (v,w) \right| \leq C_{B,1} \tbar{v} \, \tbar{w}$.

\item[(iv)]
Restricted Coercivity:
There exists constants $\gamma_0 \geq 0$ and $C_{B,2} = C_{B,2}(\gamma_0) \geq 0$,
both independent of $h$, such that, for all $\gamma \geq \gamma_0$,
there holds $B_h(v_h,v_h) \geq C_{B,2} \tbar{v_h}^2$
for all $v_h \in V_{h,r} \subset H^2(\mathcal{T}_h)$.

\item[(v)]
Equivalency:
If $\tbar{v}_B \equiv \sqrt{B_h(v,v)}$,
then
\[
	\sqrt{C_{B,2}} \tbar{v_h} \leq \tbar{v_h}_B \leq \sqrt{C_{B,1}} \tbar{v_h}
\]
for all $v_h \in V_{h,r} \subset H^2(\mathcal{T}_h)$.
\end{itemize}
\end{lemma}

\begin{remark} \rm \quad

(a)
By the Lax-Milgram theorem, the (symmetric) IPDG method
has a unique solution.
Moreover, if $u \in H^{s+1}(\Omega)$ is the solution to \eqref{BVP} with $1 \leq s \leq r$,
then there exists a positive constant $C$ independent of $h$ such that
$\| u - u_h \|_{L^2(\Omega)} \leq C h^{s+1} \| D^{s+1} u \|_{L^2(\Omega)}$.

(b)
Consistency only holds for $u \in H^2(\Omega)$.
The derivation of the formulation for $B_h$ is based on the assumption
$u \in H^2(\Omega)$ and the test functions $w \in H^2(\mathcal{T}_h)$.
The inconsistency with \eqref{BVP_weak_grad} comes from the lack of the term
$\big\langle [\nabla v] , \big\{ \varphi \big\} \mathbf{n} \big\rangle_{\mathcal{E}_h^I}$,
a term that disappears if $v \in H^2(\Omega)$.
We have also added the term
$\big\langle [v] , \{\nabla w\} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}$
in the formulation of $B_h$ where we have again used the ansatz that for $v \in H^2(\Omega)$,
the term is zero-valued.

(c)
Coercivity only holds for functions in the discrete space $V_{h,r}$.
Since the formulation is consistent, we have a solution exists in the space $H^2(\mathcal{T}_h)$.
Thus, the lack of coercivity for the entire broken Sobolev space $H^2(\mathcal{T}_h)$
implies a potential loss of uniqueness.
Notice that the lack of coercivity arises from the fact we cannot control the norm associated with
$H^2(\mathcal{T}_h)$ using the norm $\tbar{\cdot}$.
We will use distributional derivatives as a means to preserve uniqueness.

(d)
All of the results are dependent upon choosing $\gamma > 0$ sufficiently large.
In contrast, LDG methods only require the penalty constant $\gamma > 0$
to preserve existence and uniqueness.
We will explore this result further in the following sections.
\end{remark}

\section{Formal motivation}\label{heuristic_sec}

In this section, we heuristically motivate the analytic results in the following 
section. We will consider a formal definition of a ``broken" derivative operator
for broken Sobolev functions that acts as an alternative to the piecewise gradient
considered above in Sections~\ref{broken_spaces} and \ref{IPDG_sec}.
Thus, we formally let $\overline{\nabla}_h : H^1(\mathcal{T}_h) \to \mathbf{L}^2(\mathcal{T}_h)$ 
be characterized locally by
\begin{equation} \label{bdo}
	\big( \overline{\nabla}_h v , \boldsymbol{\varphi} \big)_K
	= \big\langle \{v\} , \boldsymbol{\varphi} \cdot \vec{n}_K \big\rangle_{\partial K}
		- \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_K
	\quad \forall \boldsymbol{\varphi} \in \mathbf{H}^1(K)
\end{equation}
for all $K \in \mathcal{T}_h$ for $v \in H^1(\mathcal{T}_h)$ arbitrary.
Observe, if $v \in H^1(\Omega)$, then we would have 
$\overline{\nabla}_h v = \nabla_h v = \nabla v$.

\begin{remark} \rm \quad 

(a) As seen in \cite{FengLewisNeilan13}, the operator $\overline{\nabla}_h$ forms the basis for LDG
methods when we restrict the test functions and the range of the operator to be in the finite
dimensional spaces $V_{h,r} \subset H^1(\mathcal{T}_h)$ and $\mathbf{V}_{h,r} \subset \mathbf{H}^1(\mathcal{T}_h)$,
respectively.
Thus, the results obtained in this section are a prelude to the stability results for LDG methods
that we present in Section~\ref{LDG}.

(b)
If $[ v ] \neq 0$, then $\overline{\nabla}_h v \notin \mathbf{L}^1_{\text{loc}}(\mathcal{T}_h) \supset \mathbf{L}^2(\mathcal{T}_h)$.
In general, $\overline{\nabla}_h v$ must be understood as a (potentially singular) distribution
as will be considered in the following section.
However, many of the formal results obtained below based on the assumption
$\overline{\nabla}_h v \in \mathbf{L}^2(\mathcal{T}_h)$
will be seen to carry over into the distributional setting
and provide insight into LDG methods that are based on
the restriction of \eqref{bdo} to test and trial functions in the finite dimensional space
$V_{h,r} \subset H^1(\mathcal{T}_h)$.

(c)
Since each simplex $K$ is convex and polygonal, the trace operator
$T : H^1(K) \to H^{1/2}(\partial K)$ has a bounded right inverse, \cite{grisvard}.
Thus, if $v \in H^1(\mathcal{T}_h)$ and $\{ v \} \in L^2(\partial K) \setminus H^{1/2}(\partial K)$
for some simplex $K$,
it follows that there does not exist an element $w \in H^1(\mathcal{T}_h)$ such that
$\overline{\nabla}_h v = \nabla_h w$.
If $v \in H^1(\mathcal{T}_h) \setminus C^0(\Omega)$,
then we expect such discrepancies/jumps to occur
on interior faces/edges.
Therefore, the ``broken" derivative operator $\overline{\nabla}_h$ has the
potential to address issues related to the discontinuities of a general broken Sobolev function
unlike the standard piecewise gradient.

(d)
For comparison, the piecewise gradient operator $\nabla_h : H^1(\mathcal{T}_h) \to \mathbf{L}^2(\mathcal{T}_h)$
is fully characterized locally by
\begin{equation}\label{nablah}
	\big( \nabla_h v , \boldsymbol{\varphi} \big)_K
	= \big\langle v , \boldsymbol{\varphi} \cdot \vec{n}_K \big\rangle_{\partial K}
		- \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_K
	\quad \forall \boldsymbol{\varphi} \in \mathbf{H}^1(K)
\end{equation}
for all $K \in \mathcal{T}_h$
for $v \in H^1(\mathcal{T}_h)$ arbitrary.
\end{remark}

Observe that \eqref{bdo} is equivalent to
\begin{align} \label{avg_global}
	\big( \overline{\nabla}_{h} v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
	& = - \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
		+ \big\langle \{v\} , \big[ \boldsymbol{\varphi} \big] \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
		+ \big\langle v , \boldsymbol{\varphi} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B}
\end{align}
for all $\boldsymbol{\varphi} \in \mathbf{H}^1(\mathcal{T}_h)$
while summing \eqref{nablah} over all simplices $K$ yields
\begin{align}\label{standard}
	\big( \nabla_h v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
	& = - \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
		+ \big\langle [v] , \big\{ \boldsymbol{\varphi} \big\} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
		+ \big\langle \big\{ v \big \} , \big[ \boldsymbol{\varphi} \big] \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I} \\
	\nonumber & \quad
		+ \big\langle v , \boldsymbol{\varphi} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B}
\end{align}
for all $\boldsymbol{\varphi} \in \mathbf{H}^1(\mathcal{T}_h)$.
Thus, $\overline{\nabla}_h$ is formally related to $\nabla_h$ by
\begin{align}\label{formal_compare}
	\big( \overline{\nabla}_{h} v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
	& = \big( \nabla_h v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
		- \big\langle [v] , \big\{ \boldsymbol{\varphi} \big\} 
\cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
\end{align}
for all $\boldsymbol{\varphi} \in \mathbf{H}^1(\mathcal{T}_h)$,
a relationship that will be made precise in Section~\ref{LDG}.
In contrast to the standard piecewise gradient operator $\nabla_h$,
we see that $\overline{\nabla}_h v$
implicitly incorporates jump information of the function $v$ into its definition.
The incorporation of jump information directly into the value of 
$\overline{\nabla}_h v$
is the precise reason that $\overline{\nabla}_h$ must be treated as a distribution.
The definition for $\overline{\nabla}_h$ is not entirely local;
instead, it uses the trace information of its nearest neighbors.
Therefore, we expect the ``broken" derivative operator $\overline{\nabla}_h$
to be more faithful to the global gradient operator with respect to jumps in a function,
as will be seen in the following subsections.


\subsection{Friedrichs inequality} \label{Poinc_sec}

We now derive a result analogous to Lemma~\ref{Friedrichs_standard}
that is based on the formal operator $\overline{\nabla}_h$.
We will see that we no longer require control of interior jumps
in order to maintain stability, even when considering functions in the broken
Sobolev space $H^1(\mathcal{T}_h)$.
We also note that the bound will trivially hold
in the case $\overline{\nabla}_h v \notin \mathbf{L}^2(\mathcal{T}_h)$.
Such a result is the first step in explaining the much weaker stability requirements
for LDG methods when compared to IPDG methods.

Let $v \in H^1(\mathcal{T}_h)$ and suppose $\overline{\nabla}_h v \in \mathbf{L}^2(\Omega)$.
Define $\psi \in H^2(\Omega) \cap H^1_0(\Omega)$ by $- \Delta \psi = v$.
Then, since $\Omega$ is convex,
there exists a constant $C_1$ depending only on $\Omega$ such that
$\| \psi \|_{H^2(\Omega)} \leq C_1 \| v \|_{L^2(\Omega)}$, see 
\cite{Gilbarg_Trudinger01}.
Thus, we have
\begin{align*}
\| v \|_{L^2(\Omega)}^2
& = \big( v , v \big)_{\mathcal{T}_h}
= \big( v , - \Delta \psi \big)_{\mathcal{T}_h}
= - \big( v , \nabla \cdot \nabla \psi \big)_{\mathcal{T}_h} \\
& = \big( \overline{\nabla}_h v , \nabla \psi \big)_{\mathcal{T}_h}
	- \big\langle \{v\} , \big[ \nabla \psi \big] \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
	- \big\langle v , \nabla \psi \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B} \\
& = \big( \overline{\nabla}_h v , \nabla \psi \big)_{\mathcal{T}_h}
	- \big\langle v , \nabla \psi \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B} ,
\end{align*}
where we have used the fact
$\big[ \nabla \psi \big] = 0$ in the interior of the domain.
Observe that, by Young's inequality and the trace inequality,
\begin{align*}
	\big( \overline{\nabla}_h v , \nabla \psi \big)_{\mathcal{T}_h}
	&\leq \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)} \, \| \nabla \psi \|_{L^2(\Omega)} \\
	&\leq  \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)} \, \| \psi \|_{H^2(\Omega)} \\
	&\leq C_1 \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)} \, \| v \|_{L^2(\Omega)} \\
	&\leq \frac12 C_1^2 \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)}^2
		+ \frac12 \| v \|_{L^2(\Omega)}^2
\end{align*}
and
\begin{align*}
	\big\langle v , \nabla \psi \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B}
	& \leq \sum_{e \in \mathcal{E}_h^B} \| v \|_{L^2(e)} \, \| \nabla \psi \|_{L^2(e)} \\
	& \leq \frac{\theta}{2} \sum_{e \in \mathcal{E}_h^B} h_e^{-1} \| v \|_{L^2(e)}^2
		+ \frac{1}{2 \theta} \sum_{e \in \mathcal{E}_h^B} h_e \| \nabla \psi \|_{L^2(e)}^2
\end{align*}
with
\begin{align*}
	h_e \| \nabla \psi \|_{L^2(e)}^2
	& \leq \widehat{C} \left( \| \nabla \psi \|_{L^2(K)}^2 + h_e^2 \| D^2 \psi \|_{L^2(K)}^2 \right)
	\leq C' \| \psi \|_{H^2(K)}^2 ,
\end{align*}
for $C'$ a constant such that $C' \geq \max \{\widehat{C} , \widehat{C} h_e^2 \}$
for all $e \in \mathcal{E}_h^B$.
Hence,
\begin{align*}
	\| v \|_{L^2(\Omega)}^2
	& \leq
	\frac12 C_1^2 \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)}^2
		+ \frac12 \| v \|_{L^2(\Omega)}^2
	+ \frac{\theta}{2} \sum_{e \in \mathcal{E}_h^B} h_e^{-1} \| v \|_{L^2(e)}^2  \\
	& \quad + \frac{N}{2 \theta} C' C_1 \| v \|_{L^2(\Omega)}^2
\end{align*}
for $N$ the maximum number of boundary faces/edges a boundary simplex can have.
Choosing $\theta = 2 N C' C_1$ and assuming $h_e \leq 1$ for all faces/edges, we have
\[
	\| v \|_{L^2(\Omega)}^2
	\leq C \Big( \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)}^2
		+ \sum_{e \in \mathcal{E}_h^B} h_e^{-1} \| v \|_{L^2(e)}^2 \Big) ,
\]
for $C$ independent of $h_e$ for all $e \in \mathcal{E}_h$.
Thus, we have proved the following result.


\begin{lemma} \label{Fried_cor}
There exists a constant $C$ depending only on $\Omega$ and the lower bounds for the interior
angles and adjacent-edge ratios of the triangulation such that
\[
	\| v \|_{L^2(\Omega)}
	 \leq C \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)}
\]
for all $v \in H^1(\mathcal{T}_h) \cap L^2_0(\Omega)$
and
\[
	\| v \|_{L^2(\Omega)}^2
	 \leq C \Big(  \| \overline{\nabla}_h v \|_{L^2(\mathcal{T}_h)}^2
		+ \sum_{e \in \mathcal{E}_h^B} h_e^{-1} \| v \|_{L^2(e)}^2
		\Big)
\]
for all $v \in H^1(\mathcal{T}_h)$.
\end{lemma}



\begin{remark} \rm
Observe that Lemma~\ref{Fried_cor} does not require interior jump penalization
unlike the standard result given by Lemma~\ref{poinc_standard} that is based 
on the piecewise gradient operator $\nabla_h$.
Using Theorem~\ref{stability_distribution} below, we can also conclude that
for $v \in H^1(\mathcal{T}_h)$, there holds
$\overline{\nabla}_h v = \mathbf{0}$ if and only if $v$ is a continuous, 
constant-valued function.
Thus, we would expect DG methods based on $\overline{\nabla}_h$ to not have
the same uniqueness/stability issues observed in Sections~\ref{broken_spaces} 
and \ref{IPDG_sec}.
\end{remark}

\subsection{Formal existence and uniqueness results for Poisson's equation}

We now formally extend \eqref{BVP_weak} to broken Sobolev functions
in a way that preserves existence and uniqueness properties.
The analytic extension will be given in Section~\ref{LDG} when we
reformulate \eqref{BVP_weak} using results based on distributional derivatives.

First, we define an appropriate norm so that we can apply the Lax-Milgram theorem.
Define $\tbar{\cdot}_{\overline{1}} : H^1(\mathcal{T}_h) \to [0,\infty]$ by
\[
	\tbar{v}_{\overline{1}} \equiv \big( \overline{\nabla}_h v ,
 \overline{\nabla}_h v \big)_{\mathcal{T}_h}
\]
for all $v \in H^1(\mathcal{T}_h)$.
Then, from the broken Poincar\`e/Friedrichs inequalities above, we immediately have
$\tbar{\cdot}_{\overline{1}}$ defines an (extended) norm
on $H^1(\mathcal{T}_h) \cap L^2_0(\Omega)$ and a semi-norm on $H^1(\mathcal{T}_h)$.

Using the ``broken" derivative operator $\overline{\nabla}_h$,
we formally recast \eqref{BVP_weak} as:
Find $u \in H^1(\mathcal{T}_h) \cap L^2_g(\Omega)$ such that
\begin{equation} \label{PDE_broken_form}
\big( \overline{\nabla}_h u , \overline{\nabla}_h \varphi \big)_{\mathcal{T}_h}
	= \big(f , \varphi \big)_{\mathcal{T}_h}
\end{equation}
for all $\varphi \in H^1_0(\mathcal{T}_h)$.
Since $H^1(\mathcal{T}_h)$ is a Hilbert space when equipped with the (extended) norm
$\| \cdot \|^2 \equiv \| \cdot \|_{L^2(\mathcal{T}_h)}^2 + \tbar{\cdot}_{\overline{1}}$
and coercivity follows directly from the Friedrichs' inequalities found in 
Lemma~\ref{Fried_cor},
the formal problem has a unique solution by the Lax-Milgram theorem.

\begin{remark} \quad

(a)
The above formulation does not explicitly consider jumps in $u$.
However, the formal operator $\overline{\nabla}_h$ does implicitly consider
jumps when compared to $\nabla_h$.

(b)
The above formulation is consistent with the traditional formulation in $H^1(\Omega)$
whenever the solution $u \in H^1(\Omega) \cap L^2_g(\Omega)$.
In other words, if the solution is in $H^1(\Omega)$, then it is also a formal solution to \eqref{BVP_weak}.

(c)
By using the broken derivative operator $\overline{\nabla}_h$
instead of the global gradient operator $\nabla$,
the formal problem \eqref{PDE_broken_form} posed in the broken Sobolev space $H^1(\mathcal{T}_h)$
preserves the form of the original problem \eqref{BVP_weak} posed in the 
Sobolev space $H^1(\Omega)$.
\end{remark}


\section{Distributional derivatives and local discontinuous Galerkin methods} 
\label{LDG}

We now make the observations in Section~\ref{heuristic_sec}
rigorous by extending the above results to the proper setting
based on distributional derivatives.
Distributions will allow us to extend \eqref{BVP_weak} to
the broken Sobolev space $H^1(\mathcal{T}_h)$
in such a way that consistency and uniqueness regarding the solution
$u \in H^1(\Omega) \cap L^2_g(\Omega)$
for \eqref{BVP_weak} are preserved
when testing by the larger class of functions in $H^1_0(\mathcal{T}_h)$.
The new formulation will immediately rule out any functions in
the space
$\left( H^1(\mathcal{T}_h) \cap L^2_g(\Omega) \right) \setminus 
\left( H^1(\Omega) \cap L^2(\Omega) \right)$
from being a solution, yielding a stability result that can be utilized by 
DG methods.

We first consider the distributional derivative of a broken Sobolev function.
Next we will prove a stability result that will be the key observation
motivating the use of approximate distributional derivatives
when formulating DG methods.
We will end the section by discussing the extension of our
continuous results in $H^1(\mathcal{T}_h)$ to the finite dimensional setting of DG methods.

\begin{definition} \rm
Let $v \in W^{1,p}(\mathcal{T}_h)$.
The {\it distributional derivative} of $v$ is defined as the distribution
$\mathcal{D}(v) : \mathbf{C}^\infty_c(\Omega) \to \mathbb{R}$ such that
\[
	< \mathcal{D} (v), \boldsymbol{\varphi} >
 \equiv 	- \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_\Omega
\]
for all $\boldsymbol{\varphi} \in \mathbf{C}^\infty_c(\Omega)$.
\end{definition}

\begin{remark} \rm \quad

(a)
For comparison, the distribution associated with $\nabla_h v$ is defined by
\[
	< \nabla_h (v) , \boldsymbol{\varphi} > \, \equiv
	- \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_\Omega
	+ \big\langle [v] , \boldsymbol{\varphi} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
	= \big( \nabla v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
\]
for all $\boldsymbol{\varphi} \in \mathbf{C}^\infty_c(\Omega)$.
If $v \in W^{1,p}(\Omega)$, then we have
\[
	< \nabla_h (v), \boldsymbol{\varphi} > \, = \, < \mathcal{D} (v) , \boldsymbol{\varphi} > \,
	= \, < \nabla (v) , \boldsymbol{\varphi} > \,
	\equiv - \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_\Omega
\]
for all $\boldsymbol{\varphi} \in \mathbf{C}^\infty_c(\Omega)$.
The distributional derivative
will be shown to ensure stability through the ansatz that the definition of the derivative
itself will naturally account for jumps.

(b)
The distributional derivative can be associated with a function in
$\mathbf{L}^2(\Omega)$ if and only if $v \in H^1(\Omega)$.

(c)
If we let $\boldsymbol{\varphi} \in \mathbf{C}^\infty_c(\Omega)$ in \eqref{avg_global},
then we have $\overline{\nabla}_h v$ formally satisfies
\[
	\big( \overline{\nabla}_h v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
	= - \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_{\Omega}
\]
for all $v \in W^{1,p}(\mathcal{T}_h)$.
Thus, the formal operator considered in Section~\ref{heuristic_sec}
behaves like a regularized/idealized version of the distribution $\mathcal{D}$.
\end{remark}

\subsection{A stability result for distributional derivatives}

A property of functions in $H^1(\Omega)$ is that
$\| \nabla v \|_{L^2(\Omega)} = 0$ if and only if $v$ is constant valued 
for all $v \in H^1(\Omega)$.
We now extend the result to distributional derivatives over the larger
function space $H^1(\mathcal{T}_h)$.
As noted above, such a result does not hold when considering $\nabla_h$.
Thus, $\nabla_h$ has a larger nullspace than $\mathcal{D}$ yielding a loss 
of uniqueness
when deriving a weak formulation of \eqref{BVP}.

\begin{theorem} \label{stability_distribution}
Let $v \in H^1(\mathcal{T}_h)$.
Then $v$ is a continuous, constant-valued function over $\Omega$
if and only if
$\mathcal{D}(v) = \vec{0}$,
where $\vec{0}$ is understood as the zero-valued function in the set of all
bounded, linear functionals acting on $\mathbf{C}^\infty_c(\Omega)$.
\end{theorem}

\begin{proof}
The forward direction is immediate.
Let $w$ be a constant valued function over $\Omega$.
Since $w \in H^1(\Omega)$, we have $\mathcal{D}(w) = \nabla w = \vec{0}$.

We now prove the reverse direction.
We use a variational proof that is based on choosing appropriate test functions
since such a proof is commonplace in the DG literature.
A much shorter proof is given in Remark~\ref{distribution_remark}
that uses results from distribution theory.

Let $v \in H^1(\mathcal{T}_h)$ and suppose $\mathcal{D} (v) = \vec{0}$.
We first show $v$ is piecewise constant with respect to the triangulation.
To this end, we will choose suitable test functions that ensure $\nabla_h v = \mathbf{0}$.
Pick $K \in \mathcal{T}_h$.
Let $\eta_{K,\epsilon} \in C^\infty(\overline{K})$ denote a ``cutoff'' function on $K$,
\cite{Gilbarg_Trudinger01},
such that
\[
	\eta_{K,\epsilon}(x) = \begin{cases}
	1 & \text{if } \operatorname{dist}(x,\partial K) \geq 2 \, \epsilon , \\
	0 & \text{if } \operatorname{dist}(x,\partial K) \leq \epsilon
	\end{cases}
\]
for $\epsilon > 0$ sufficiently small.
Let $\rho_\epsilon \in C^\infty \big( B(0,\epsilon) \big)$ denote the standard
mollifier function on $B(0,\epsilon)$, see \cite{Evans}.
Then, we have $\operatorname{supp}(\rho_\epsilon) \subset B(0,\epsilon)$.
Define $K_\epsilon = \{ x \in K \mid \operatorname{dist} (x,\partial K) \geq \epsilon \}$.
Then, $\rho_\epsilon * \varphi \in C^\infty( K_\epsilon )$ for all $\varphi \in L^2(\Omega)$,
where $*$ denotes the convolution operator.
Notationally, we let $\varphi^\epsilon \equiv \rho_\epsilon * \varphi$.
Analogous results hold for $\boldsymbol{\varphi} \in \mathbf{L}^2(\Omega)$
and $\boldsymbol{\varphi}^\epsilon \equiv \rho_\epsilon \mathbf{1} * \boldsymbol{\varphi}$.

By the definition of $\mathcal{D} (v)$, we have
\[
	0 =  < \mathcal{D}(v) , \boldsymbol{\varphi}  >
	= - \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_{\Omega}
\]
for all $\boldsymbol{\varphi} \in \mathbf{C}^\infty_c(\Omega)$.
Letting $\boldsymbol{\varphi} = \eta_{K,\epsilon} 
 \nabla v^\epsilon \in \mathbf{C}^\infty_c(K)$,
where the $\epsilon$ superscript denotes a mollified version of the function,
we have
$\operatorname{supp} ( \boldsymbol{\varphi} ) \subset K_\epsilon$ and
\begin{align*}
	0 & = - \big( v , \nabla \cdot (\eta_{K,\epsilon} \, \nabla v^\epsilon) \big)_{K} \\
	& = \big( \nabla v , \eta_{K,\epsilon} \, \nabla v^\epsilon \big)_K
		- \big\langle v , \eta_{K,\epsilon} \, \nabla v^\epsilon \cdot \vec{n}_K \big\rangle_{\partial K} \\
	& = \big( \nabla v , \eta_{K,\epsilon} \, \nabla v^\epsilon \big)_K \\
	& = \big( \nabla v , \nabla v^\epsilon \big)_{K_{2 \epsilon}}
		+ \big( \nabla v , \eta_{K,\epsilon} \, \nabla v^\epsilon
			 \big)_{K_\epsilon \setminus K_{2 \epsilon}} .
\end{align*}
Since $\nabla v^\epsilon \to \nabla v$ in $L^2(K_\epsilon)$, $\nabla v \in L^2(K)$, and
$\mu \left( K_\epsilon \setminus K_{2 \epsilon} \right) \to 0$,
we have $\nabla v \big|_K = \mathbf{0}$.
Therefore, $v$ is constant on $K$.

We now show $v \in C^0(\Omega)$.
Pick $e \in \mathcal{E}_h^I$.
Let $K^+, K^- \in \mathcal{T}_h$ such that $\partial K^+ \cap \partial K^- = e$,
and let $x^+, x^-$ denote the barycentric centers of $K^+, K^-$, respectively.
Define $\Omega_e \subset \mathbb{R}^d$ as the convex domain corresponding to
the simplex formed by the vertices of $e$ and the nodes $x^+$ and $x^-$.

By the definition of $\mathcal{D} (v)$ and the fact that $v$ is piecewise constant
with respect to the triangulation, there holds
\begin{align*}
	0 & = - \big( v , \nabla \cdot \boldsymbol{\varphi} \big)_{\Omega}
	= \big( \nabla_h v , \boldsymbol{\varphi} \big)_{\mathcal{T}_h}
		- \big\langle [v] , \boldsymbol{\varphi} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
	= - \big\langle [v] , \boldsymbol{\varphi} \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
\end{align*}
for all $\boldsymbol{\varphi} \in C^\infty_c(\Omega)$.
Let $j \in \{1,2,\ldots,d \}$ and let
$\boldsymbol{\varphi} = \varphi \, \mathbf{e}_j$
for $\varphi \in C^\infty_c(\Omega)$ and $\mathbf{e}_j$ the $j$th Cartesian basis vector.
Then,
\[
	0 = - \big\langle [v] , \varphi \, n_j \big\rangle_{\mathcal{E}_h^I} .
\]

Let $\eta_{e,\epsilon} \in C^\infty(\overline{e})$
denote a ``cutoff'' function on $e$ such that
\[
	\eta_{e,\epsilon}(x) = \begin{cases}
	1 & \text{if } \operatorname{dist}(x,\partial e) \geq 2 \, \epsilon , \\
	0 & \text{if } \operatorname{dist}(x,\partial e) \leq \epsilon
	\end{cases}
\]
for $\epsilon > 0$ sufficiently small.
Define $e_{\epsilon} = \{ x \in e \mid \operatorname{dist} (x,\partial e) \geq \epsilon \}$.
Let $\rho'_\epsilon \in C^\infty \big( B(0,\epsilon) \big)$ denote the standard
mollifier function on $B(0,\epsilon)$ in $\mathbb{R}^{d-1}$.
Choose $\varphi \in C^\infty_c(\Omega)$ such that
$\varphi \big|_e = \eta_{e , \epsilon} \, \rho'_\epsilon * [v] \, \text{sign } n_j$ and
$\operatorname{supp} \varphi \subset \Omega_e$,
where we have used the fact that $\mathbf{n}$ is piecewise constant on $\mathcal{E}_h^I$
and $\overline{e} \cap \mathring{e}_1 = \emptyset$ (in the topology for $\mathbb{R}^{d-1}$)
for all $e_1 \in \mathcal{E}_h$ such that $e \neq e_1$.
Then, we have
\[
	\big\langle [v] ,
		\eta_{e , \epsilon} \, \rho'_\epsilon * [v] \, | n_j | \big\rangle_e
	= 0
\]
for all $j = 1, 2, \ldots, d$.
Letting $\epsilon \to 0$ and using the fact $[v] \in L^2(e)$,
we have $v$ is continuous across $e$.
Therefore, $v$ is a continuous, constant-valued function over $\Omega$
since $v$ is both piecewise constant and continuous.
\end{proof}

\begin{remark} \label{distribution_remark} \rm
Using results from distribution theory, we could identify the zero distribution uniquely
(up to sets of measure zero)
as the zero function $\mathbf{0} \in \mathbf{L}^1_{\text{\rm loc}}(\Omega)$.
Since $\mathbf{0} \in \mathbf{L}^2(\Omega)$,
we have $v \in H^1(\Omega)$ with $\nabla v = \mathbf{0}$.
Thus, $v$ is a continuous, constant-valued function over $\Omega$.
Note that the underlying results for distributions in this short argument
are typically proved using density arguments combined with mollifiers,
as was done in the proof of the theorem.
\end{remark}

We can now precisely extend \eqref{BVP_weak} to broken Sobolev spaces
by exploiting properties of distributional derivatives.
We will have that the solution $u$ to problem \eqref{BVP_weak}
is also a solution to the new generalized problem.
Furthermore, the solution to the generalized problem will be unique due to
Theorem~\ref{stability_distribution}.
Thus, the new generalized problem posed in broken Sobolev spaces
extends the existence and uniqueness properties of the weak formulation
given by \eqref{BVP_weak}.

In general, the distribution $\mathcal{D}(v)$ associated with $v \in W^{1,p}(\mathcal{T}_h)$
cannot be associated with a function $\mathcal{D}v \in \mathbf{L}^1_{\text{loc}}(\Omega)$ defined by
\[
	\big( \mathcal{D} v , \boldsymbol{\varphi} \big)_\Omega
	= < \mathcal{D}(v) , \boldsymbol{\varphi} >
\]
for all $\boldsymbol{\varphi} \in \mathbf{C}^\infty_c(\Omega)$.
The problem is that for discontinuous functions in $W^{1,p}(\mathcal{T}_h)$,
the distributional derivative acts like a delta function/measure
which has no representation in $\mathbf{L}^1_{\text{loc}}(\Omega)$, \cite{Evans}.
Such a negative result is the key to successfully extending
\eqref{BVP_weak} to broken Sobolev spaces without a loss of uniqueness.

Let $\varphi \in H^1_0(\mathcal{T}_h)$ and
$\left\{ \varphi_k \right\}_{k=1}^\infty \subset C^\infty_c(\Omega)$
such that $\varphi_k$ converges weakly to $\varphi$ with respect to the
$L^2$ inner product, i.e., $\varphi_k \to^* \varphi$
(cf. \cite[Theorem 6.32]{Rudin}).
Then,
if $u \in H^1(\Omega) \cap L^2_g(\Omega)$, we have
\begin{align*}
	\big( f , \varphi_k \big)_\Omega
	& = \big(\mathcal{D} u , \mathcal{D} \varphi_k \big)_\Omega
	= \big( \nabla u , \mathcal{D} \varphi_k \big)_\Omega
	= \big( \nabla u , \nabla \varphi_k \big)_\Omega
\end{align*}
for all $k \geq 1$.
Using the fact $\varphi_k \to^* \varphi$, we can conclude
\[
	\big( f , \varphi \big)_\Omega = \lim_{k \to \infty} \big( \nabla u , \nabla \varphi_k \big)_\Omega ,
\]
consistent with the weak formulation \eqref{BVP_weak}.

Without a loss of generality, we can assume $g = 0$ and $f \neq 0$.
Then, if $u \in H^1_0(\mathcal{T}_h) \setminus H^1(\Omega)$, we trivially have there exists a simplex $K$
such that
\[
	\big( f , u \big)_K
	\neq \big(\mathcal{D} u , \mathcal{D} u \big)_K ,
\]
where we have let $\varphi = u \chi_K$ for $\chi_K$
the characteristic function supported on $K$.
Thus, the problem:
Find $u \in H^1(\mathcal{T}_h) \cap L^2_g(\Omega)$ such that
\begin{equation} \label{bvp_dist}
	\big(\mathcal{D} u , \mathcal{D} \varphi \big)_\Omega
	= \big( f , \varphi \big)_\Omega
\end{equation}
for all $\varphi \in H^1_0(\mathcal{T}_h)$
successfully captures the essentials of \eqref{BVP_weak}
without introducing new solutions in the larger discontinuous space $H^1(\mathcal{T}_h) \cap L^2_g(\Omega)$,
where the left-hand side of the formulation must be interpreted in an approximate sense
based on density arguments.


\subsection{Application to local discontinuous Galerkin methods}

We first introduce a simple version of the LDG method
that will form as the basis for the application of results from 
Sections~\ref{heuristic_sec} and \ref{LDG}.
To this end, we introduce the central DG derivative operator that will be used
in the formulation of the LDG method.
Define the spaces
\[
	C^m(\mathcal{T}_h) \equiv \prod_{K \in \mathcal{T}_h} C^m \left( \, \overline{K} \, \right),
	\quad \quad
	\mathcal{V}_h \equiv W^{1,1}(\mathcal{T}_h)\cap C^0(\mathcal{T}_h)
\]
for $m \geq 0$ an integer.
Notice $V_{h,r} \subset \mathcal{V}_h$.

\begin{definition} \label{bdoh} \rm
The \emph{central DG derivative operator}
$\overline{\nabla}_{h,r} : \mathcal{V}_h \to \mathbf{V}_{h,r}$ is defined locally by
\begin{equation*}
	\big( \overline{\nabla}_{h,r} v , \boldsymbol{\varphi}_h \big)_K
	\equiv \big\langle \{v\} , \boldsymbol{\varphi}_h \cdot \vec{n}_K \big\rangle_{\partial K}
		- \big( v , \nabla \cdot \boldsymbol{\varphi_h} \big)_K
	\quad \forall \boldsymbol{\varphi}_h \in \mathbf{V}_{h,r}
\end{equation*}
for all $K \in \mathcal{T}_h$ for all $v \in \mathcal{V}_h$.
Then, summing over all simplexes $K \in \mathcal{T}_h$, we have
\[
	\big( \overline{\nabla}_{h,r} v , \boldsymbol{\varphi}_h \big)_{\mathcal{T}_h}
	= \big\langle \{v\} , [\boldsymbol{\varphi}_h] \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
		+ \big\langle v , \boldsymbol{\varphi}_h \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B}
		- \big( v , \nabla \cdot \boldsymbol{\varphi_h} \big)_{\mathcal{T}_h}
\]
for all $\boldsymbol{\varphi}_h \in \mathbf{V}_{h,r}$.
\end{definition}


Notice that every function $v \in \mathcal{V}_h$
has a well-defined trace in $L^1(\partial K)$
and every function $\boldsymbol{\varphi}_h \in \mathbf{L}^\infty(\partial K)$ 
for all $K \in \mathcal{T}_h$. Thus, by the auxiliary result in 
\cite[Appendix A]{LewisNeilan12}, we have
the central DG derivative operator is well-defined
for all functions in $\mathcal{V}_h$.
If $v \in H^1(\Omega)$, then we immediately have the following approximation result.

\begin{theorem}[\cite{FengLewisNeilan13}]\label{prop0a}
For any $v\in \mathcal{V}_h\cap H^1(\Omega)$,
$\overline{\nabla}_h v$ coincides with the $L^2$ projection of $\nabla v$
onto $\mathbf{V}_{h,r}$.
We write $\overline{\nabla}_h v = \mathcal{\vec{P}}_{h,r} \nabla v$,
where $\mathcal{\vec{P}}_{h,r}$ denotes the $L^2$ projection 
operator onto $\mathbf{V}_{h,r}$.
\end{theorem}

\begin{remark} \rm \quad

(a)
Formally, we have $\overline{\nabla}_{h,r}$ is the $L^2$ projection of the
operator $\overline{\nabla}_h$ considered in Section~\ref{heuristic_sec}.

(b)
For $v_h \in V_{h,r}$, we have $\overline{\nabla}_{h,r} v_h \in \mathbf{V}_{h,r-1}$
if $v_h$ is also in $H^1(\Omega)$.
Furthermore, there exist functions $v_h \in V_{h,r}$ such that
$\overline{\nabla}_{h,r} v_h \notin \mathbf{V}_{h,r-1}$.
For comparison, we have the piecewise gradient operator satisfies
$\nabla_h \left( V_{h,r} \right) \subset \mathbf{V}_{h,r-1}$.

(c)
If the function $v \in \mathcal{V}_h$ has given Dirichlet boundary data,
we define the \emph{central DG derivative operator with given Dirichlet boundary data}
$\overline{\nabla}^g_{h,r} : \mathcal{V}_h \to \mathbf{V}_{h,r}$ by
\[
	\big( \overline{\nabla}^g_{h,r} v , \boldsymbol{\varphi}_h \big)_{\mathcal{T}_h}
	= \big\langle \{v\} , [\boldsymbol{\varphi}_h]
\cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^I}
		+ \big\langle g , \boldsymbol{\varphi}_h \cdot \mathbf{n} \big\rangle_{\mathcal{E}_h^B}
		- \big( v , \nabla \cdot \boldsymbol{\varphi_h} \big)_{\mathcal{T}_h}
\]
for all $\boldsymbol{\varphi}_h \in \mathbf{V}_{h,r}$.
\end{remark}

Finally, using the results found in \cite{FengLewisNeilan13},
we can define the LDG method (in primal form) by:
find $u_h \in V_{h,r}$ such that
\begin{equation} \label{ldg_primal}
	\big( \overline{\nabla}^g_{h,r} u_h , \overline{\nabla}^{0}_{h,r} \varphi_h \big)_{\mathcal{T}_h}
	+ \gamma \big\langle [u_h] , \big[ \varphi_h \big] \big\rangle_{\mathcal{E}_h^I}
	+ \gamma \big\langle u_h - g , \varphi_h \big\rangle_{\mathcal{E}_h^B}
	= \big( f , \varphi_h \big)_{\mathcal{T}_h}
\end{equation}
for all $\varphi_h \in V_{h,r}$.
Then, for $r \geq 1$ and $\gamma > 0$, the LDG method defined by \eqref{ldg_primal}
has a unique solution (see \cite{Cockburn98,unified}).
Thus, the proposed method based on $\overline{\nabla}_{h,r}$ is stable as long as any
positive penalty parameter is used.
Unlike the (symmetric) IPDG method defined in Section~\ref{IPDG_sec} based on the
piecewise gradient operator $\nabla_h$,
we do not need to require the penalty parameter be sufficiently large in order to
guarantee uniqueness.
If we further assume $\gamma = \alpha h^{-1}$ for an appropriate positive 
constant $\alpha$
chosen independently of $h$
and $u \in H^{r+2}(\Omega)$ for $u$ the solution to \eqref{BVP_weak},
then there holds
\[
	\| u - u_h \|_{L^2(\Omega)}
	+ h \| \nabla u - \overline{\nabla}^g_{h,r} u_h \|_{L^2(\Omega)}
	\leq C h^{r+1} \| u \|_{H^{r+2}(\Omega)}
\]
for some positive constant $C$ independent of $h$ (cf. \cite{Castilo00}).

\begin{remark} \rm
If we choose $X_{h,r} \subset V_{h,r} \cap H^1(\Omega)$
and add the restriction $u_h , \varphi_h \in X_{h,r}$,
then \eqref{ldg_primal} is equivalent to the finite element method
with weakly enforced boundary conditions.
\end{remark}

The formal stability results found in Section~\ref{heuristic_sec}
hint that only boundary penalization should be sufficient
for the LDG method.
However, as seen in Figure~\ref{LDG_example},
for particular choices of $r$ and $h$
it is possible to construct non-constant, discontinuous functions
in the nullspace of $\overline{\nabla}^0_{h,r}$
yielding a loss of uniqueness for the LDG method when choosing $\gamma = 0$.
Thus, the LDG method given by \eqref{ldg_primal} does require $\gamma > 0$
due to the loss of information when projecting onto the finite dimensional space.
However, the MD-LDG method developed by Cockburn and Dong in \cite{CockburnDong07}
has extended the need for penalization to only a subset of the boundary faces/edges,
consistent with the Friedrichs inequality given by Lemma~\ref{Fried_cor}.
The main idea is to choose appropriate interior trace values in Definition~\ref{bdoh}
that eliminate the need to penalize along interior faces/edges.
All of the interior trace values are in the set $\big\{\{ v \} , v^+ , v^- \big\}$.
By incorporating two complementing DG derivative operators into the
discretization of \eqref{BVP} and averaging,
nonstandard LDG methods were developed in \cite{LewisNeilan12},
called dual-wind DG methods,
that eliminated the need for penalization along both interior and boundary 
faces/edges.
Thus, the boundary conditions appear naturally in the method
based on the choice of the underlying DG derivative operators.
A short survey of when various DG methods do not require interior 
and/or boundary penalization
can be found in Section 2.2 of \cite{feng_lewis_wise15}.

\begin{figure} \label{LDG_example}
\begin{center}

\begin{tikzpicture}[scale = 2]
\draw (0,0) rectangle (3,3);

\draw (0,3) -- (3,0);
\draw (0,2) -- (2,0);
\draw (0,1) -- (1,0);
\draw (1,3) -- (3,1);
\draw (2,3) -- (3,2);
\draw (1,0) -- (1,3);
\draw (2,0) -- (2,3);
\draw (0,1) -- (3,1);
\draw (0,2) -- (3,2);

\node at (0.1,0.1) {1};
\node at (0.1,0.75) {-1};
\node at (0.25,0.9) {1};
\node at (0.1,1.1) {-1};
\node at (0.1,1.75) {0};
\node at (0.25,1.9) {0};
\node at (0.1,2.1) {0};
\node at (0.1,2.75) {1};
\node at (0.25,2.9) {-1};

\node at (1.1,0.1) {0};
\node at (1.1,0.75) {1};
\node at (1.25,0.9) {-1};
\node at (1.1,1.1) {1};
\node at (1.1,1.75) {-1};
\node at (1.25,1.9) {1};
\node at (1.1,2.1) {-1};
\node at (1.1,2.75) {0};
\node at (1.25,2.9) {0};

\node at (2.1,0.1) {-1};
\node at (2.1,0.75) {0};
\node at (2.25,0.9) {0};
\node at (2.1,1.1) {0};
\node at (2.1,1.75) {1};
\node at (2.25,1.9) {-1};
\node at (2.1,2.1) {1};
\node at (2.1,2.75) {-1};
\node at (2.25,2.9) {1};

\node at (0.75,0.1) {0};
\node at (0.9,0.25) {0};
\node at (0.9,0.9) {-1};
\node at (0.75,1.1) {1};
\node at (0.9,1.25) {-1};
\node at (0.9,1.9) {1};
\node at (0.75,2.1) {-1};
\node at (0.9,2.25) {1};
\node at (0.9,2.9) {0};

\node at (1.75,0.1) {-1};
\node at (1.9,0.25) {1};
\node at (1.9,0.9) {0};
\node at (1.75,1.1) {0};
\node at (1.9,1.25) {0};
\node at (1.9,1.9) {-1};
\node at (1.75,2.1) {1};
\node at (1.9,2.25) {-1};
\node at (1.9,2.9) {1};

\node at (2.75,0.1) {1};
\node at (2.9,0.25) {-1};
\node at (2.9,0.9) {1};
\node at (2.75,1.1) {-1};
\node at (2.9,1.25) {1};
\node at (2.9,1.9) {0};
\node at (2.75,2.1) {0};
\node at (2.9,2.25) {0};
\node at (2.9,2.9) {-1};

\end{tikzpicture}

\end{center}
\caption{
Example of a function in the nullspace of $\overline{\nabla}^0_{h,1}$, 
cf. \cite{brezzi}.
The function is piecewise linear with the value of $1$, $-1$, or $0$ at 
each node in the mesh.
Note that the function is not in the nullspace of $\overline{\nabla}_{h,1}$.}
\end{figure}


We end this section by relating the LDG method to our results concerning 
distributional derivatives.
We show that $\overline{\nabla}_{h,r}$ is in fact the $L^2$ projection of
$\mathcal{D}(v)$ onto the finite dimensional space $\mathbf{V}_{h,r}$
for all $v \in \mathcal{V}_h$,
when interpreted correctly.
Proofs of the following results can all be found in \cite{FengLewisNeilan13}.

Let $v \in \mathcal{V}_h$, $\mathcal{D} (v)$ denote the distributional
derivative of $v$,
and $\Xi \subset \Omega$ be a ($d$-1)-dimensional continuous and bounded surface.
We define the \emph{delta function} $\delta(\Xi, g, x)$
 of variable strength supported on $\Xi$ by
(cf. \cite{Engquist_Tornberg_Tsai05})
\begin{align}\label{e4.0a}
	< \delta (\Xi, g, x), \varphi >
	\equiv \int_\Xi g(s) \varphi(x(s))\, ds
\end{align}
for all $\varphi\in C^0(\Omega)$,
where $x(s) \in \Xi$.
We also extend the above definition to
test functions from $V_{h,r}$ by
\begin{align}\label{e4.0b}
	< \delta (\Xi, g, x), \varphi_h > \,
	\equiv \int_\Xi g(s) \big\{\varphi_h(x(s))\big\} \, ds
\end{align}
for all $\varphi_h\in V_{h,r}$.
Using $\delta(\Xi, g, x)$, we have $\mathcal{D}(v)$ is characterized by
\begin{align}\label{e4.1b}
	\mathcal{D} (v)
	= \sum _{K \in \mathcal{T}_h} \nabla v \, \chi_K
		- \sum_{e \in \mathcal{E}_h^I} \mathbf{n} \, \delta(e, [v], x) 
\quad \mbox{for a.e. } x \in \Omega,
\end{align}
where $\chi_K$ denotes the characteristic function supported on $K$.
Thus, we are able to express the central DG derivative operator
as the projection of the distributional derivative of $v$ as follows:

\begin{theorem}\label{prop0b}
For any $v \in \mathcal{V}_h$,
$\overline{\nabla}_{h,r} v$ coincides with the $L^2$ projection of
$\mathcal{D} (v)$ onto $\mathbf{V}_{h,r}$ in the sense that
\begin{align}\label{e4.1e}
	\big( \overline{\nabla}_{h,r} v, \boldsymbol{\varphi}_h \big)_{\mathcal{T}_h}
		= \, < \mathcal{D} (v), \boldsymbol{\varphi}_h >
\end{align}
for all $\boldsymbol{\varphi}_h \in \mathbf{V}_{h,r}$,
where the right hand-side is understood according to \eqref{e4.0b}.
We write $\overline{\nabla}_{h,r} v = \mathcal{\vec{P}}_{h,r} \mathcal{D} (v)$.
\end{theorem}

Using the above theorem, we can also formally compare
the central DG derivative operator to the piecewise gradient operator
using a lifting operator as follows:

\begin{corollary}
Define the ``lifting operator" $\mathcal{L}_{h,r}: L^1(\mathcal{E}_h) \to \mathbf{V}_{h,r}$ by
\begin{equation}\label{lift}
	\big( \mathcal{L}_{h,r} v, \boldsymbol{\varphi}_h \big)_{\mathcal{T}_h}
	\equiv \sum_{e \in \mathcal{E}_h^I} 
< \mathbf{n}(e)  \delta(e, [v], \cdot),  \boldsymbol{\varphi}_h >
	= \big \langle [v] \mathbf{n}, \big\{ \boldsymbol{\varphi}_h \big\} \big\rangle_{\mathcal{E}_h^I}
\end{equation}
for all $\boldsymbol{\varphi}_h \in \mathbf{V}_{h,r}$.
Then we have
\begin{equation}\label{lift_1}
	\overline{\nabla}_{h,r} v = \nabla_h v - \mathcal{L}_{h,r} v
\end{equation}
for all $v \in L^1(\mathcal{E}_h)$.
\end{corollary}

Thus, we see that jump stabilization is inherently enforced by the lifting operator
that arises when comparing $\overline{\nabla}_{h,r}$ to the piecewise
gradient operator $\nabla_h$.


Finally, combining Theorems \ref{prop0a} and \ref{prop0b}
and the well-known limiting characterization theorem of distributional
derivatives (cf. \cite[Theorem 6.32]{Rudin}), we obtain another
characterization for our central DG derivative operator.

\begin{theorem}\label{prop0c}
For any $v \in \mathcal{V}_h$,
there exists a sequence of functions
$\{v_j \}_{j\geq 1} \subset C^\infty_c(\Omega)$
such that
\begin{enumerate}
\item[(i)]
$v_j \to  v$ as $j \to \infty$ in $L^1(\Omega)$.
\item[(ii)]
$\overline{\nabla}_{h,r} v_j \to \mathcal{\vec{P}}_{h,r} \mathcal{D} (v)$
as $j \to \infty$ weakly in $\mathbf{V}_{h,r}$ in the sense that
\begin{equation}\label{e4.1g}
	\lim_{j \to \infty} \big( \overline{\nabla}_{h,r} v_j, \boldsymbol{\varphi}_h \big)_{\mathcal{T}_h}
	= \big(\mathcal{\vec{P}}_{h,r} \mathcal{D} (v), \boldsymbol{\varphi}_h \big)_{\mathcal{T}_h}
\end{equation}
for all $\boldsymbol{\varphi}_h \in \mathbf{V}_{h,r}$.
\end{enumerate}
\end{theorem}


Therefore, we have the LDG method is equivalent to a finite dimensional
restriction of the distributional formulation of \eqref{BVP_weak} over
 broken Sobolev spaces given by \eqref{bvp_dist},
where the Dirichlet boundary data is enforced weakly
and interior penalization is added to ensure the preservation of uniqueness
when projecting onto the finite dimensional solution space.

We make one more observation concerning the relationship between
distributional derivatives and flux-based DG methods written in primal form.
By adjusting the extension of the delta function of variable strength for
test functions in $V_{h,r}$ given by \eqref{e4.0b},
we are able to recover other DG derivative operators.
The adjustment is simply choosing different flux values for $\varphi_h$
such as $\varphi_h^+$ or $\varphi_h^-$.
Thus, we can extend all of the above results and interpret a wide class of DG methods
(including the MD-LDG method)
as simply the $L^2$ projection of the broken Sobolev formulation \eqref{bvp_dist}
onto finite dimensional solution spaces
where penalization is added on interior and/or boundary faces/edges
only when necessary to ensure uniqueness and/or the enforcement of boundary conditions.
We also see that DG methods that better approximate the distributional derivative
instead of the piecewise gradient operator inherit better stability properties
explaining the fundamental difference between the (symmetric) IPDG method 
and the LDG method with regards to stabilization.


\section{Conclusion} \label{conclusion_sec}

The underlying goal of this paper was to extend uniqueness results for second 
order elliptic PDEs to broken Sobolev spaces
as motivation for the design and understanding of various DG approximation methods.
We showed that IPDG methods are based on a formulation that does not preserve
uniqueness when the problem is posed in a broken Sobolev space
using the natural choice of a piecewise gradient operator.
In contrast, LDG methods are based on a formulation that guarantees uniqueness
by considering distributional derivatives,
derivative operators that naturally incorporate jump information.
By combining uniqueness and consistency results for distributional derivatives,
we see that not only does a solution exist in the broken Sobolev space setting,
but the solution corresponds to the solution of the original PDE problem posed 
in (full) Sobolev spaces.
Therefore, using distributional derivatives to extend PDE problems to broken
 Sobolev spaces provides the correct context for motivating a wide class of 
DG finite element approximation methods based on projections.

Another consequence of the above results is
that flux-based DG methods can actually be considered
conformal methods for the PDE problem \eqref{BVP}
when the PDE problem is understood in terms of distributional derivatives.
By posing the weak formulation \eqref{BVP_weak} in the larger space of 
broken Sobolev functions,
both continuous Galerkin methods (such as the finite element method)
and DG methods can be discussed with respect to the same unifying PDE setting.
The choice of (discrete) solution space $X_{h,r} \subset V_{h,r}$
and the choice of the DG derivative operator fully determine the 
approximation method and whether or not penalization
is necessary for a given mesh $\mathcal{T}_h$.
Therefore, the results found in this paper naturally
motivate the use of the weak Galerkin finite element method (\cite{wang2013weak}) 
and the DG finite element differential calculus (\cite{FengLewisNeilan13})
for understanding and designing stable DG methods for approximating solutions
to partial differential equations.

\begin{thebibliography}{00}

\bibitem{unified} D. Arnold, F. Brezzi, B. Cockburn,  D. Marini;
\emph{Unified analysis of discontinuous {G}alerkin methods for elliptic problems},
SIAM J. Numer. Anal., 39:1749--1779, 2001.

\bibitem{Baker77}  G. A. Baker;
\emph{Finite element methods for elliptic equations using nonconforming elements},
Math. Comp., 31:45--59, 1977.

\bibitem{bassi_rebay}  F. Bassi, S.Rebay;
\emph{A high-order accurate discontinuous finite element method for the numerical 
solution of the compressible Navier-Stokes equations},
J. Comput. Phys., 131:267--279, 1997.

\bibitem{brenner03}  S. C. Brenner;
\emph{Poincar\`e-Friedrichs inequalities for piecewise $H^1$ functions},
SIAM J. Numer. Anal., 41(1): 306--324, 2003.

\bibitem{Brenner}  S. C.~Brenner and L.R.~Scott;
\emph{The Mathematical Theory of Finite Element Methods} (Third edition),
Springer, 2008.

\bibitem{brezzi}  F.~Brezzi, M.~Manzini, D.~Marini, P.~Pietra, A.~Russo;
\emph{Discontinuous finite elements for Diffusion problems},
Francesco Brioschi (1824-1897) Convegno di Studi Matematici,
October 22-23, 1997,
Ist. Lomb. Acc. Sc. Lett., Incontro di studio N. 16, 197-217, 1999.

\bibitem{Castilo00}  P.~Castillo, B.~Cockburn, I.~Perugia, D.~Sch\"otzau;
\emph{An a priori error analysis of the local discontinuous Galerkin method 
for elliptic problems}, SIAM J. Numer. Anal., 38(5): 1676--1706, 2000.

\bibitem{Ciarlet78}  P.~G.~Ciarlet;
\emph{The Finite Element Method for Elliptic Problems},
North-Holland, Amsterdam, 1978.

\bibitem{CockburnDong07}  B.~Cockburn and B.~Dong;
\emph{An analysis of the minimal dissipation local discontinuous
Galerkin method for convection-diffusion problems},
J. Sci. Comput., 32(2): 233--262, 2007.

\bibitem{cockburn2009unified}
 B. Cockburn,  J. Gopalakrishnan,  R. Lazarov;
\emph{Unified hybridization of discontinuous Galerkin, mixed, and continuous
 Galerkin methods for second order elliptic problems},
SIAM J. Numer. Anal., 47: 1319--1365, 2009.

\bibitem{Cockburn98}  B.~Cockburn, C-W.~Shu;
\emph{The local discontinuous Galerkin
method for time-dependent convection-diffusion systems},
SIAM J. Numer. Anal., 35(6): 2440--2463, 1998.

\bibitem{DouglasDupont76}  J. Douglas Jr., T. Dupont;
\emph{Interior Penalty Procedures for Elliptic and Parabolic
Galerkin Methods}, Lecture Notes in Phys. 58, Springer-Verlag, Berlin, 1976.

\bibitem{Engquist_Tornberg_Tsai05}
 B. Engquist, A. K. Tornberg,  R. Tsai;
\emph{Discretization of dirac delta functions in level set methods},
J Comput. Phy., 207:28--51, 2005.

\bibitem{Evans}  L.~C. Evans;
\emph{Partial Differential Equations, Graduate Studies in Mathematics}, 2nd Edition,
American Mathematical Society, Providence, RI, 2010.

\bibitem{feng_lewis_wise15}  W. Feng, T. Lewis, and S. Wise;
\emph{Discontinuous Galerkin derivative operators with applications to
second order elliptic problems and stability},
 Mathematical Meth. in App. Sciences, DOI: 10.1002/mma.3440, 2015.

\bibitem{FengLewisNeilan13}  X.~Feng, T.~Lewis, M.~Neilan;
\emph{Discontinuous Galerkin
finite element differential calculus and applications to numerical solutions
of linear and nonlinear partial differential equations},
accepted by J. Comput. Appl. Math.,
 arXiv:1302.6984 [math.NA].

\bibitem{Gilbarg_Trudinger01}  D.~Gilbarg, N.~S. Trudinger;
\emph{Elliptic Partial Differential Equations of Second Order,
 Classics in Mathematics},
Springer-Verlag, Berlin, 2001, reprint of the 1998 edition.

\bibitem{grisvard}  P.~Grisvard;
\emph{Elliptic Problems in Nonsmooth Domains},
Pitman, Boston, MA, 1985.

\bibitem{LewisNeilan12}  T.~Lewis, M.~Neilan;
\emph{Convergence analysis of a symmetric dual-wind discontinuous Galerkin method},
J. Sci. Comput., 59(3):602--625, 2014.

\bibitem{RiviereWheelerGirault00}  B.~Rivi\`ere, M.F.~Wheeler, V.~Girault;
\emph{Improved energy estimates for interior penalty,
constrained and discontinuous Galerkin methods for elliptic problems. I.},
 Comput. Geosci., 3(3--4): 337--360, 1999.

 \bibitem{RiviereWheelerGirault01}
  B.~Rivi\`ere, M. F.~Wheeler,  V.~Girault;
\emph{A priori error estimates for finite element methods based
on discontinuous approximation spaces for elliptic problems},
SIAM J. Numer. Anal. 39(3): 902--931, 2001.

\bibitem{Rudin}  W. Rudin;
\emph{Functional Analysis}, Second Edition, McGraw-Hill, New York, 1991.

\bibitem{wang2013weak}  J. Wang, X. Ye;
\emph{A weak Galerkin finite element method for second-order elliptic problems},
J. Comput. Appl. Math., 241:103--115, 2013.

\bibitem{Wheeler78}  M. F. Wheeler;
\emph{An elliptic collocation-finite element method with interior penalties},
SIAM J. Numer. Anal., 15: 152--161, 1978.

\end{thebibliography}

\end{document}