path: root/stresstensor.tex

\chapter{Construction of the Stress Tensor of~a~Free~Scalar~Quantum~Field}
\label{chapter:stress-tensor}

\begin{center}
  \emph{Note: Work in Progress}
\end{center}

\begin{equation*}
  H = \tfrac{1}{2} \parens*{(\partial_t \phi)^2 + \abs{\nabla_{\!\!\symbfit{x}} \phi}^2 + m^2 \phi^2} 
\end{equation*}

At the end of this \namecref{chapter:stress-tensor} 
we will have gained the ability to rigorously define
arbitrary renormalized products of the free field and its derivatives
as a densely defined quadratic-form valued tempered distribution,
which on the dense subspace of the smooth vectors of the Hamiltonian
is realized by essentially selfadjoint operators.

\section{Choosing Conventions and Fixing Notation}
\label{section:conventions}

It is an unfortunate reality of quantum physics literature that
there is a great deal of variation in notation and choice of signs and constants.
While this does not affect the physical or mathematical content,
it is a hindrance when working with formulas from multiple sources.
In the present \namecref{section:conventions} we detail our choices
as a service to the reader.

\begin{itemize}
  \item \emph{Minkowski space} $M=\RR^4$ equipped with the \emph{Lorentz bilinear form} (or metric)
    \begin{equation*}
      x \cdot y = g_{\mu \nu} x^{\mu} y^{\nu} = x^0y^0 - x^1 y^1 - x^2 y^2 - x^3 y^3
    \end{equation*}
    points $x = (x^0,x^1,x^2,x^3) \in M$ are sometimes written $x = (x^0,\symbfit{x})$ with separated time and space coordinates
  \item Given a complex-valued function $f$ on $M$,
    we define its \emph{Fourier transform}\index{Fourier transform} $\ft{f}\,$ by
    \nomenclature[f]{$\ft{f}$}{Fourier transform of $f$}
    \begin{equation}
      \label{fourier-transform}
      \ft{f}(p) \defequal \int_{M} e^{i p \cdot x} f(x) \, dx
    \end{equation}
    whenever the integral converges. The \emph{inverse Fourier transform} is defined by
    \begin{equation*}
      \label{inverse-fourier-transform}
      \ift{f}(p) \defequal \frac{1}{(2 \pi)^2} \int_{M} e^{-i p \cdot x} f(x) \, dx.
    \end{equation*}
  \item To a mathematician $\overline{\phantom{z}}$ usually means complex conjugation and ${}^*$ indicates the Hilbert adjoint of an operator,
    while a physicist may read ${}^*$ as complex conjugation and
    denotes the Hilbert adjoint with ${}^{\dagger}$.
    We avoid confusion by using $\overline{\phantom{z}}$ for complex conjugation
    and ${}^{\dagger}$ for the Hilbert adjoint exclusively.
  \item $\schwartz{M}$ is the space of complex-valued Schwartz functions on $M$ \\
    $\realschwartz{M}$ is the space of real-valued Schwartz functions on $M$ 
  \item The \emph{Bosonic Fock space} over a Hilbert space $\hilb{H}$ is denoted $\BosonFock{\hilb{H}}$. \\
    Its \emph{finite particle subspace} is denoted $\BosonFockFinite{\hilb{H}}$.
  \item Abstract free field: The \emph{Segal quantization} $\Phi$ assigns to every
    $g \in \hilb{H}$, a selfadjoint (unbounded) operator $\Phi(g)$ in $\BosonFock{\hilb{H} }$,
    which on the the finite particle subspace is given by
    \begin{equation*}
      \Phi_{\mathrm{S}} (g) = \frac{1}{\sqrt{2}} \parens*{a(g) + a(g)^\dagger} 
    \end{equation*}
    annihilation and creation operators, $g \in \hilb{H}$, $\psi \in \BosonFock{\hilb{H}}$ for $\hilb{H} = L^2(R^4,\Omega_m)$
    \begin{align*}
      \parens[\big]{a(g) \psi} {}_n (k_1, \ldots, k_n)
      &= \sqrt{n+1} \int_M \! \bar{g}(p) \, \psi_{n+1} (p,k_1, \ldots, k_n) \, d\Omega_m(p) \\
      \parens[\big]{a(g)^\dagger \psi} {}_n (k_1, \ldots, k_n)
      &= \frac{1}{\sqrt{n}} \sum_{i=1}^n g(k_i) \, \psi_{n-1} (k_1, \ldots, \widehat{k_i}, \ldots, k_n)
    \end{align*}
    The symbol $\widehat{\hphantom{k_i}}$ over $k_i$ indicates omission.
    \begin{align*}
      \parens[\big]{a(g) \psi} {}_{n-1} (k_1, \ldots, k_{n-1}) 
      &= \sqrt{n} \int_M \! \bar{g}(p) \, \psi_n (p,k_1, \ldots, k_{n-1}) \, d\Omega_m(p) \\
      \parens[\big]{a(g)^\dagger \psi} {}_{n+1} (k_1, \ldots, k_{n+1}) 
      &= \frac{1}{\sqrt{n+1}} \sum_{i=1}^{n+1} g(k_i) \, \psi_n (k_1, \ldots, \widehat{k_i}, \ldots, k_{n+1})
    \end{align*}
    \begin{multline*}
      \parens[\big]{a(g) a(g) \psi} {}_{n-2} (k_1, \ldots, k_{n-2}) = \\
      \sqrt{n} \sqrt{n-1} \int_M \int_M \! \bar{g}(p_1) \bar{g}(p_2) \, \psi_n (p_1,p_2,k_1, \ldots, k_{n-s}) \, d\Omega_m(p_1) d\Omega_m(p_2) \\
    \end{multline*}
    For later use, we also give the action of an $s$-fold product of annihilation operators:
    \begin{multline}
      \label{equation:multiple-annihilation-operators}
      \parens[\big]{a(g_1) \cdots a(g_s) \psi} {}_{n-s} (k_1, \ldots, k_{n-s}) = 
      \sqrt{n (n-1) \cdots (n-s+1)} \cdot {} \\
      \cdot \int_M \!\! d\Omega_m(p_1) \cdots \!\! \int_M \!\! d\Omega_m(p_s) \ \bar{g_1}(p_1) \cdots \bar{g_s}(p_s) \ \psi_n (k_1, \ldots, k_{n-2},p_1,\ldots,p_s)
    \end{multline}
  \item We consider the free Hermitian scalar field of mass $m > 0$. \\
    \emph{mass hyperboloid} $X_m^+ = \braces{p \in M \mid p^2 = m^2, p^0 > 0 }$ 
    with normalized Lorentz invariant measure $\Omega_m$
  \item single particle state space: $\hilb{H} = L^2(X_m^+, \Omega_m)$
    \begin{equation*}
      E : \schwartz{M} \to \hilb{H}, \quad f \mapsto Ef = \left.\ft{f}\,\right\vert {X_m^+}
    \end{equation*}
    We define a $\RR$-linear mapping $\phi$ by
    \begin{equation*}
      \realschwartz{M} \ni f \mapsto \varphi(f) = \Phi_{\mathrm{S}}(Ef) = \frac{1}{\sqrt{2}} \parens*{a(Ef) + a(Ef)^\dagger}
    \end{equation*}
    This extends to complex valued test functions $f \in \schwartz{M}$
    \begin{equation*}
      \varphi(f) = \frac{1}{\sqrt{2}} \parens*{a(E\bar{f}) + a(Ef)^\dagger}
    \end{equation*}
    called the \emph{massive free scalar quantum field} 
  \item
    annihilation and creation operators, $f \in \schwartz{M}$, $\psi \in \BosonFock{\hilb{H}}$
    \begin{align*}
      \parens[\big]{a(f) \psi} {}_n (k_1, \ldots, k_n)
      &= \sqrt{n+1} \int_M \! \overline{Ef(p)} \, \psi_{n+1} (p,k_1, \ldots, k_n) \, d\Omega_m(p) \\
      \parens[\big]{a(f)^\dagger \psi} {}_n (k_1, \ldots, k_n)
      &= \frac{1}{\sqrt{n}} \sum_{i=1}^n Ef(k_i) \, \psi_{n-1} (k_1, \ldots, \widehat{k_i}, \ldots, k_n)
    \end{align*}
\end{itemize}

\begin{equation*}
  \normord{\varphi(f)^2} = \tfrac{1}{2} \parens[\big]{a^{\dagger}(Ef) a^{\dagger}(Ef) + a(Ef) a(Ef)} + a^{\dagger}(Ef) a(Ef) 
\end{equation*}

\begin{equation*}
  \normord{\varphi(f)^2} = \tfrac{1}{2} \parens[\big]{a^{\dagger} a^{\dagger} + a a} + a^{\dagger} a \quad \text{where} \quad a = a(Ef)
\end{equation*}

\begin{equation*}
  \innerp{\psi'}{\normord{\varphi(f)^2} \,\psi} = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \innerp{\psi'_m}{\normord{\varphi(f)^2} \,\psi_n}
\end{equation*}

There will be no contribution to the sum
unless either $m=n+2$ or $m=n$ or $m=n-2$.
Now let us calculate those contributions.
To avoid the introduction of sums,
we use the adjoint identity
to transform creation operators on the right
into annihilation operators on the left.

\begin{align*}
  \innerp{\psi'_{n+2}}{\normord{\varphi(f)^2} \,\psi_n} 
  &= \tfrac{1}{2} \innerp{a(Ef) a(Ef) \psi'_{n+2}}{\psi_n} \\
  &= \begin{multlined}[t][10cm]
    \tfrac{1}{2} \sqrt{n+2} \sqrt{n+1} \int dp_1 dp_2 \, \ft{f}(p_1) \ft{f}(p_2) \int dk_1 \cdots dk_n \\
    \overline{\psi'_{n+2}(p_1,p_2,k_1,\ldots,k_n)} \, \psi_n(k_1,\ldots,k_n)
  \end{multlined} \\[1ex]
  \innerp{\psi'_{n}}{\normord{\varphi(f)^2} \,\psi_n} 
  &= \innerp{a(Ef) \psi'_{n}}{a(Ef) \psi_n} \\
  &= \begin{multlined}[t][10cm]
    \sqrt{n} \sqrt{n} \int dp_1 dp_2 \, \ft{f}(p_1) \overline{\ft{f}(p_2)} \int dk_1 \cdots dk_{n-1} \\
    \overline{\psi'_{n}(p_1,k_1,\ldots,k_{n-1})} \, \psi_n(p_2,k_1,\ldots,k_{n-1})
  \end{multlined} \\[1ex]
  \innerp{\psi'_{n-2}}{\normord{\varphi(f)^2} \,\psi_n} 
  &= \tfrac{1}{2} \innerp{\psi'_{n-2}}{a(Ef) a(Ef) \psi_n} \\
  &= \begin{multlined}[t][10cm]
    \tfrac{1}{2} \sqrt{n} \sqrt{n-1} \int dp_1 dp_2 \, \overline{\ft{f}(p_1)} \overline{\ft{f}(p_2)} \int dk_1 \cdots dk_{n-2} \\
    \overline{\psi'_{n-2}(k_1,\ldots,k_{n-2})} \psi_n(p_1,p_2,k_1,\ldots,k_{n-2})
  \end{multlined}
\end{align*}

\begin{align*}
  \innerp{\psi'_{m}}{\normord{\varphi(f)^2} \,\psi_n} 
  &= \begin{multlined}[t][10cm]
    \tfrac{1}{2} \sqrt{n+2} \sqrt{n+1} \int dp_1 dp_2 \, \ft{f}(p_1) \ft{f}(p_2) \int dk_1 \cdots dk_n \\
    \chi(p_{s+1}) \cdots \chi(p_2) \\
    \overline{\psi'_{m}\parens{k_1,\ldots,k_{n-s},p_1,\ldots,p_s}} \, \psi_n\parens{k_1,\ldots,k_{n-(2-s)},p_{s+1},\ldots,p_2}
  \end{multlined}
\end{align*}

\section{Quadratic Forms}
\label{section:quadratic-forms}

In a typical physics literature treatment of second quantization,
the annihilation and creation operators and the quantized field are treated as
point-dependent operator valued functions $a(p)$, $a^\dagger(p)$, $\varphi(x)$,
disregarding the fact that these may not be operators, in a strict sense, and without smearing with a test function.
Nonetheless, this notational fiction is useful, and we can uphold it with little effort by giving
the pointwise \enquote{operators} rigorous meaning as quadratic forms.

Let us consider a real scalar quantum field with mass parameter $m>0$.
Thus, the single-particle state space is the Hilbert space $\hilb{H} = L^2(X_m^+,\Omega_m)$,
where $X_m^+$ is the upper half of the hyperboloid defined by the condition $p \cdot p = m^2$ in momentum space,
and $\Omega_m$ is the unique normalized Lorentz-invariant measure on it.
Given a point $p$ in momentum space,
we define the annihilation operator $a(p)$ to be the operator
whose domain is the finite particle subspace $\BosonFockFinite{\hilb{H}}$ of Bosonic Fock space
and whose action on a vector $\psi = (\psi_n)_{n \ge 0}$ in $\BosonFockFinite{\hilb{H}}$ is given by
\begin{equation*}
  \parens[\big]{a(p) \psi} {}_n (k_1, \ldots, k_n)
  = \sqrt{n+1} \, \psi_{n+1} (k_1, \ldots, k_n,p).
\end{equation*}
While this is a perfectly well-defined operator,
an issue arises when one looks for an adjoint to this operator.
A formal calculation based on the adjoint identity
\begin{equation}
  \label{equation:adjoint-identity}
  \innerp[\big]{\psi'}{a(p)^\dagger \psi} =
  \innerp[\big]{a(p) \psi'}{\psi}
\end{equation}
leads to
\begin{equation}
  \label{equation:creation-operator-at-point}
  \parens[\big]{a(p)^\dagger \psi} {}_n (k_1, \ldots, k_n)
  = \frac{1}{\sqrt{n}} \sum_{i=1}^n \delta(p - k_i) \, \psi_{n-1} (k_1, \ldots, \widehat{k_i}, \ldots, k_n),
\end{equation}
where the symmetrization is necessary to obtain an expression that
at least has a chance of being a $n$ Boson state.
However, it clearly is not a $L^2$ function.
Given any state $\psi'$, we can
formally calculate the inner product of $\psi'$ with~\eqref{equation:creation-operator-at-point} 
and use the result to define $a^\dagger(p)$
as a mapping that assigns a number to each \emph{pair} of states.
That is, we define the creation \enquote{operator} $a^\dagger(p)$
to be the quadratic form
\begin{gather*}
  a(p)^\dagger \vcentcolon \BosonFockFinite{\hilb{H}} \times \BosonFockFinite{\hilb{H}} \longrightarrow \CC \\
  \innerp[\big]{\psi'}{a(p)^\dagger \psi}
  \defequal
  \begin{multlined}[t]
    \sum_{n=1}^{\infty} \frac{1}{\sqrt{n}} \sum_{i=1}^n
    \int_M \! d\Omega_m(k_1) \cdots \widehat{d\Omega_m(k_i)} \cdots d\Omega_m(k_n) \\
    \cdot \overline{\psi'_{n} (k_1, \ldots, \underset{i}{p}, \ldots, k_n)}
    \psi_{n-1} (k_1, \ldots, \widehat{k_i}, \ldots, k_n).
  \end{multlined}
\end{gather*}
One can verify directly that
with this definition the adjoint identity~\eqref{equation:adjoint-identity} 
holds for all $\psi,\psi' \in \BosonFockFinite{\hilb{H}}$.
For completeness, we give a precise definition of quadratic form.

\begin{definition}{Quadratic Form}{quadratic-form}
  A \emph{quadratic form}\index{quadratic form} $q$ on a complex Hilbert space $\hilb{H}$ is a mapping
  \begin{equation*}
    q \vcentcolon D(q) \times D(q) \to \CC,
  \end{equation*}
  where $D(q)$ is a linear subspace of $\hilb{H}$, called the \emph{form domain}\index{form domain}\index{quadratic form!domain of a},
  such that $q$ is conjugate linear in its first argument
  and linear in its second argument (i.e.\ sesquilinear).
  We say that $q$ is \emph{densely defined}
  if $D(q)$ is dense in $\hilb{H}$.
\end{definition}

The creation \enquote{operator} $a^\dagger(p)$ we considered above
is a densely defined quadratic form on Bosonic Fock space
in the sense of \cref{definition:quadratic-form}
and its form domain is the finite particle subspace.

Any linear operator on a complex Hilbert space $\hilb{H}$ has an
obvious interpretation as a quadratic form on $\hilb{H}$,
and the form domain agrees with the domain of the operator.
The reverse construction is always possible for densely defined quadratic forms,
but one may obtain an operator with trivial domain.

\begin{definition}{Operator Associated to a Quadratic Form}{}
  Suppose $q$ is a densely defined quadratic form on a complex Hilbert space $\hilb{H}$.
  The linear \emph{operator associated to}\index{quadratic form!operator associated to a} $q$, denoted $\QFop{q}$,
  is defined on the domain
  \begin{equation*}
    D(\QFop{q}) = \braces{\psi \in D(q) \mid \text{the map $q(\cdot,\psi) \vcentcolon D(q) \to \CC$ is bounded}},
  \end{equation*}
  and maps $\psi \in D(\QFop{q})$ to the vector $\QFop{q}\psi$ in $\hilb{H}$ satisfying
  $q(\psi'\!,\psi) = \innerp{\psi'\!}{\QFop{q}\psi}$,
  which exists and is unique by Riesz’s Representation Theorem.
\end{definition}

We will use the symbol $\QFequal$ between quadratic forms or operators
to indicate their equality as quadratic forms.
\todo{statement about domains?}

A natural question is how the smeared operators relate to the pointwise ones.
The answer is simple:

\begin{equation*}
  a(g) \QFequal \int \overline{g(p)} a(p) \, d\Omega_m(p)
\end{equation*}

\begin{equation*}
  a^\dagger(g) \QFequal \int g(p) a^\dagger(p) \, d\Omega_m(p)
\end{equation*}
We have to explain what is meant by the integral on the right hand side.
Suppose $q(p)$ is a quadratic form on $\BosonFock{\hilb{H}}$ for each $p \in \RR^4$,
that share a common domain $D \subset D(q(p))$,
and $g$ is in $\hilb{H} = L^2(\RR^4,\Omega_m)$. Then we define
a quadratic form by
\begin{equation*}
  \parens[\bigg]{\int g(p) q(p) \, d\Omega_m(p)}(\psi',\psi)
  \defequal \int g(p) \parens[\big]{q(p)}(\psi',\psi) \, d \Omega_m(p)
\end{equation*}
for all $\psi,\psi' \in D$.

A product of the form
\begin{equation*}
  a(p_1)^\dagger \cdots a(p_s)^\dagger a(p_{s+1}) \cdots a(p_r)
\end{equation*}
can be rigorously defined as quadratic form by setting
\begin{equation*}
  \innerp[\big]{\psi'}{a(p_1)^\dagger \cdots a(p_s)^\dagger a(p_{s+1}) \cdots a(p_r) \psi}
  \defequal \innerp[\big]{a(p_1) \cdots a(p_s) \psi'}{a(p_{s+1}) \cdots a(p_r) \psi}.
\end{equation*}

\section{Normal Ordering}
% The Renormalization Map?

\blockcquote{Wick1950}{%
 \textelp{} we then proceed to rearrange such a product so as to carry all
 creation operators to the left of all destruction operators \textelp{}. The
 main problem to be solved in carrying out this idea is one of algebraic
 technique \textelp{} 
}
The process of renormalizing a product of field operators
has the purpose of discarding infinite constants
that occur when calculating the vacuum expectation value.

In physics texts the standard prescription for doing this goes as follows:
(1) Express all field operators by creation and annihilation operators
and expand the expression into a sum of products.
(2) Within each product, move all creation operators to the left of all annihilation operators
while ignoring their commutation relations.

This ad-hoc way of dealing with infinities is not entirely satisfactory.
For once, it leaves the uniqueness question (whether there is another way to
achieve the desired effect) open.
Moreover, it obscures the facts that we are renormalizing in relation to the vacuum state
(in principle, one could use another state for reference), and that the process
is representation-independent and purely algebraic.

With this in mind, let us extract the algebraic essence of the situation.
The objects of our calculations are the field operators $\Phi(f)$,
but it does not matter that these are realized as linear maps on Fock space.
Forming the product $\Phi(f)\Phi(g)$ might as well be done purely symbolically,
since none of what we want to do depends on this product
having the meaning of operator composition;
and similar for the other two arithmetic operations,
addition and multiplication with a complex scalar.
Thus we should calculate with abstract objects $\Phi(f)$ labeled by Hilbert space vectors $f \in \hilb{H}$.
Considering that here the symbol $\Phi$ carries no meaning, we can drop it and use the label $f$ itself to represent the object.
This leads us to consider formal expressions
\begin{equation*}
  \alpha^{(0)} e + \sum_{i} \alpha^{(1)}_i z^{(1)}_i + \sum_{j,k} \alpha^{(2)}_{j,k} \, z^{(2)}_j z^{(2)}_k + \cdots
\end{equation*}
where the $z^{(1)}_i,z^{(2)}_j,z^{(2)}_k,\ldots$ are in $\hilb{H}$,
the $\alpha^{(0)},\alpha^{(1)}_i,\alpha^{(2)}_{j,k},\ldots$ are complex numbers,
of which only finitely many are nonzero,
and $e$ is a special object representing an empty product of $z$'s.
To make this mathematically precise:
we are speaking of the unital non-commutative associative algebra over $\CC$
freely generated by the elements of $\hilb{H}$.
The unit of the algebra is $e$.

This is not quite what we want,
as we yet need to account for the commutation relations
$\bracks{\Phi(f),\Phi(g)} = i \Imag \innerp{f}{g}$.
By abstract algebra, this is viable
by forming the quotient of the free algebra
with respect to the two-sided ideal
generated by all elements $zz' - z'z = i \Imag \innerp{z}{z'} \cdot e$,
where $z,z' \in \hilb{H}$.
In addition, we must implement the $\RR$-linear dependence of $\Phi(f)$ on $f$.

\begin{definition}{Infinitesimal Weyl Algebra}{}
  Let $\hilb{H}$ be a complex Hilbert space.
  The \emph{infinitesimal Weyl algebra}\index{infinitesimal Weyl algebra} $\InfinitesimalWeylAlg(\hilb{H})$ over $\hilb{H}$
  is the unital non-commutative associative algebra over $\CC$
  generated by the elements of $\hilb{H}$, with the relations
  \begin{align*}
    \alpha \cdot z &= 1 \cdot (\alpha z) \qquad \alpha \in \RR, z \in \hilb{H}, \\
    zz' - z'z &= i \Imag \innerp{z}{z'} \, e \qquad z,z' \in \hilb{H},
  \end{align*}
  where $e$ is the unit of the algebra.
  Moreover, the mapping $z \mapsto z^* = z$
  extends to a $*$-operation that turns $\InfinitesimalWeylAlg(\hilb{H})$ into a $*$-algebra.
\end{definition}
This construction is an instance of what is known as a quantization functor.
For a longer discussion of its functorial properties we refer the reader to~\cite{Fewster2012}.

To avoid any confusion, we emphasize that 
$1 \cdot iz$ (coefficient $1$, vector $iz$) and $i \cdot z$ (coefficient $i$, vector $z$)
are always different algebra elements (unless $z=0$), and there no rule that allows one to move imaginary numbers from the coefficient part to the vector part.

Now let $\Phi = \Phi_{\mathrm{S}}$ be
the Segal quantization for the free Boson field over $\hilb{H}$,
which was discussed in the preceding section.
For each $z \in \hilb{H}$,
we consider $\Phi(z)$ as an everywhere-defined linear operator
on the finite particle space $\BosonFockFinite{\hilb{H}}$ by restriction.
The mapping $z \mapsto \Phi(z)$ thus induces a natural $*$-representation
of the infinitesimal Weyl algebra,
\begin{align*}
  \Phi \vcentcolon \InfinitesimalWeylAlg(\hilb{H}) &\longrightarrow \End \parens[\big]{\BosonFockFinite{\hilb{H}}} \\
  w = z_1 \!\cdots z_n &\longmapsto \Phi(w) = \Phi(z_1) \circ \cdots \circ \Phi(z_n),
\end{align*}
since $\InfinitesimalWeylAlg(\hilb{H})$ is generated by elements of $\hilb{H}$
and $\bracks{\Phi(z),\Phi(z')} = i \Imag \innerp{z}{z'}$ on $\BosonFockFinite{\hilb{H}}$ for each $z,z' \in \hilb{H}$.

In general, any renormalization scheme for the infinitesimal Weyl algebra $\InfinitesimalWeylAlg$
will be in relation to a given linear functional, i.e.\ an element $E$ of the dual space $\InfinitesimalWeylAlg'$,
which we interpret to be yielding expectation values.
Each element of $\InfinitesimalWeylAlg$ should then be mapped to an element of $\InfinitesimalWeylAlg$
that is renormalized in the sense that it has zero expectation value in relation to $E$.
We will see that this can be accomplished by bringing each product into a special order called normal order.
In the case of the free Boson field,
the Fock vacuum state $\FockVacuum$ gives rise to the linear functional
\begin{equation*}
  \InfinitesimalWeylAlg(\hilb{H}) \ni w \mapsto E(w) \defequal \innerp[\big]{\FockVacuum}{\Phi(w) \FockVacuum},
\end{equation*}
which we call the \emph{normal vacuum}.
The normal vacuum has two essential properties: Firstly,
\begin{equation}
  \label{equation:normal-vaccum-1}
  E(e) = 1
\end{equation}
since $\Phi(e)$ is the identity operator and the vacuum $\Omega$ is a unit vector.
For stating the second property we need to introduce some notation:
\begin{definition}{Annihilator and Creator}{}
  Suppose $\InfinitesimalWeylAlg(\hilb{H})$ is the infinitesimal Weyl algebra
  over some complex Hilbert space $\hilb{H}$.
  For all $z \in \hilb{H}$
  we define, as elements of $\InfinitesimalWeylAlg(\hilb{H})$, the \emph{annihilator}
  \begin{equation*}
    \weylannihilator(z) = \frac{1}{\sqrt{2}} \parens{z+i \cdot iz},
  \end{equation*}
  and the \emph{creator}
  \begin{equation*}
    \weylcreator(z) = \frac{1}{\sqrt{2}} \parens{z-i \cdot iz}.
  \end{equation*}
\end{definition}
The second important property of the normal vacuum $E$ is
\begin{equation}
  \label{equation:normal-vaccum-2}
  E\parens[\Big]{\prod_{i=1\vphantom{S}}^{s} \weylcreator(z_i)
  \prod_{\mathclap{j=s+1\vphantom{S}}}^{r} \weylannihilator(z_j)} = 0
  \qquad \forall z_1,\ldots,z_r \in \hilb{H}, r \ge 1, 1 \le s \le r.
\end{equation}
This follows from $\Phi(\weylannihilator(z)) = a(z)$, $\Phi(\weylcreator(z)) = a^{\dagger}(z)$.
%XXX more detail

It is easy to verify that for all $z,z' \in \hilb{H}$
\begin{gather*}
  z = \frac{1}{\sqrt{2}} \parens[\big]{\weylannihilator(z) + \weylcreator(z)}, \\
  \bracks[\big]{\weylannihilator(z),\weylcreator(z')} = i \Imag \innerp{z}{z'} e.
\end{gather*}
Taken together, these identities show that $\InfinitesimalWeylAlg$ is linearly generated by expressions of the form
$\prod_{i=1}^{s} \weylcreator(z_i) \prod_{j=s+1}^{r} \weylannihilator(z_j)$.
We conclude that the normal vacuum is the \emph{unique} functional on $\InfinitesimalWeylAlg$
satisfying the conditions~\eqref{equation:normal-vaccum-1} and~\eqref{equation:normal-vaccum-2}, which are physically justified.

A \emph{monomial} in the Weyl algebra $\InfinitesimalWeylAlg$ over a complex Hilbert space $\hilb{H}$ is an element of the form
$z_1 \!\cdots z_r \in \InfinitesimalWeylAlg$, where $r \ge 0$ and $z_1,\ldots,z_r$ are in $\hilb{H}$.
We allow $r=0$, meaning that the unit $e$ is a monomial.
The set of all monomials in $\InfinitesimalWeylAlg$ is denoted $\Mon(\InfinitesimalWeylAlg)$.

\begin{definition}{Normal Ordering}{}
  Let $\hilb{H}$ be a complex Hilbert space and $\InfinitesimalWeylAlg$ its associated infinitesimal Weyl algebra.
  The mapping $\normord{\,\,}$\nomenclature[:]{$\normord{\,\,}$}{normal ordering},
  defined by linear extension of the mapping
  \begin{gather}
    \Mon(\InfinitesimalWeylAlg) \longrightarrow \InfinitesimalWeylAlg \nonumber\\
    \label{equation:normal-ordering}
    \normord{z_1 \!\cdots z_r} =
    \frac{1}{\sqrt{2^r}}
    \sum_{I \subset \braces{1,\ldots,r}} \,
    \prod_{i \in I\vphantom{\lbrace\rbrace}} \weylcreator(z_i)
    \prod_{\mathclap{j \in \braces{1,\ldots,r} \setminus I}} \weylannihilator(z_j),
  \end{gather}
  is called the \emph{normal} (or \emph{Wick}) \emph{ordering}\index{normal ordering}\index{Wick ordering} on $\hilb{H}$.
  An element $w \in \InfinitesimalWeylAlg$ is said to be in \emph{normal} (or \emph{Wick}) \emph{order},
  if $\normord{w} = w$.
\end{definition}

The products in~\eqref{equation:normal-ordering} are well defined
because creators commute with creators and annihilators commute with annihilators.
Since an empty product equals, per convention, the neutral element of multiplication, which here is the unit $e$,
the formula makes sense even for $r=0$ and asserts that $\normord{e} = e$.
The cases $r=1$ and $r=2$ read
\begin{align*}
  \normord{z} &=
  \frac{1}{\sqrt{2}} \parens[\big]{\weylannihilator(z) + \weylcreator(z)} = z, \\
  \normord{z_1 z_2} &= \frac{1}{2}
  \parens[\big]{\weylannihilator(z_1) \weylannihilator(z_2) + \weylcreator(z_1) \weylannihilator(z_2)
  + \weylcreator(z_2) \weylannihilator(z_1) + \weylcreator(z_1) \weylcreator(z_2)}.
\end{align*}
This suggests that the normally ordered product $\normord{z_1 \!\cdots z_r}$
is symmetric in $z_1,\ldots,z_n$. This is in fact true, and becomes evident
if one brings~\eqref{equation:normal-ordering} into the equivalent form
\begin{gather}
  \label{equation:normal-ordering-symmetric}
  \normord{z_1 \!\cdots z_r} =
  \frac{1}{\sqrt{2^r}}
  \sum_{\sigma \in \SymmetricGroup{r}}
  \sum_{s=0\vphantom{S}}^{r}
  \frac{1}{s!(r-s)!}
  \prod_{i=1\vphantom{S}}^{s} \weylcreator(z_{\sigma(i)})
  \prod_{\mathclap{j=s+1\vphantom{S}}}^{r} \weylannihilator(z_{\sigma(j)}) 
\end{gather}
by basic combinatorial arguments:
\nomenclature[S]{$\SymmetricGroup{r}$}{symmetric group on the set $\Set{1,\ldots,r}$}
It is easiest to start from~\eqref{equation:normal-ordering-symmetric}
and observe that the product remains invariant
when we permute creators with creators and annihilators with annihilators.
There are precisely $s!(r-s)!$ such permutations.
Let $\sigma$ be any of these.
The associated product corresponds to that summand of~\eqref{equation:normal-ordering}
for which $I = \Set{\sigma(1),\ldots,\sigma(s)}$.
The term $1 / s!(r-s)!$ cancels out.
In more technical language, we are using the fact that the mapping
\begin{align*}
  \SymmetricGroup{r} \times \Set{1,\ldots,r} &\longrightarrow \PowerSet[\big]{\Set{1,\ldots,r}} \\
  (\sigma,s) &\longmapsto \Set{\sigma(1),\ldots,\sigma(s)}
\end{align*}
is surjective and that the preimage of each $s$-element set has cardinality $s!(r-s)!$.
\nomenclature[PX]{$\PowerSet{X}$}{power set of a set $X$}
In~\cite{Klein1973}, the factor $1 / s!(r-s)!$ is erroneously missing
(which makes no difference for $r < 2$).

As a direct consequence of~\eqref{equation:normal-vaccum-2} we obtain
\begin{equation*}
  E(\normord{z_1 \!\cdots z_r}) = 0 \qquad \forall z_1,\ldots,z_r \in \hilb{H}, r \ge 1
\end{equation*}
as desired.

The normal ordered product is supposed to represent the identical quantity as before ordering,
except that we have adjusted our point of reference, so that the vacuum expectation value is zero.
It is therefore \emph{physically reasonable} to demand that the commutation relations
of the normal ordered product with the field remain unchanged.
As it turns out, this additional property makes the construction of normal ordering
\emph{mathematically unique}.
This is the content of the following theorem.

\begin{theorem}{Characterization of the Normal Order}{}
  Let $\hilb{H}$ be a complex Hilbert space and
  let $E$ be the normal vacuum on $\InfinitesimalWeylAlg(\hilb{H})$.
  Then, normal ordering is the unique mapping $N : \Mon(\InfinitesimalWeylAlg) \to \InfinitesimalWeylAlg$ such that
  \begin{gather*}
    E\parens[\big]{N(z_1 \!\cdots z_r)} = 0 \\
    \bracks{N(z_1 \!\cdots z_r), z'} =
    \sum_{s=1}^{r} \bracks{z_s,z'} N(z_1 \!\cdots \widehat{z_i} \cdots z_r)
  \end{gather*}
  for all $z_1,\ldots,z_r,z' \in \hilb{H}$ and all $r \ge 1$.
\end{theorem}
A proof of this statement is contained in the proof of~\cite[Theorem 7.1]{Baez1992}.

%\begin{theorem}{}{}
  %The normal ordering is the renormalization with respect to the normal vacuum.
%\end{theorem}

\section{Linear Differential Operators and their Formal Adjoint}

Before we turn to the problem of defining renormalized products of a quantum field and its derivatives
we must clarify what is meant mathematically by the derivative of a field.
For this, we recall that in Wightman’s approach to Quantum Field Theory,
a quantum field $\varphi$ on a spacetime manifold $M$ is modeled by an operator valued tempered distribution,
that is, a mapping that assigns to each (Schwartz class) test function $f$ on $M$ an unbounded operator $\varphi(f)$
in the Fock space $\BosonFock{\hilb{H}}$ over some Hilbert space $\hilb{H}$, such that for each fixed pair of states $\psi,\psi'$
the mapping
\begin{equation*}
  \schwartz{M} \to \CC, \quad
  f \mapsto \innerp{\psi'}{\varphi(f) \psi}
\end{equation*}
is a (scalar-valued) tempered distribution on $M$.
It is well known that tempered distributions have partial derivatives of any order.
Suppose we work with $M = \RR^d$ for simplicity,
and let $\partial_i$ denote the partial derivative with respect to the $i$-th coordinate.
Then a general \emph{linear differential operator with constant coefficients} on $M$ looks like
\begin{equation*}
  D = \sum_{\alpha} a_{\alpha} \partial^{\alpha},
\end{equation*}
where the sum runs over all multi-indices $\alpha = (\alpha_1,\ldots,\alpha_d) \in \NN^d$,
the coefficients $a_{\alpha}$ are complex numbers,
and $\partial^{\alpha} = \partial_1^{\alpha_1} \!\cdots \partial_d^{\alpha_d}$.
Then the \emph{distributional derivative} of a tempered distribution $\eta \in \tempdistrib{\RR^d}$
is defined by
\begin{equation*}
  (D\eta)(f) = \eta(D^{\dagger}f) \quad \forall f \in \mathcal{S},
\end{equation*}
where the \emph{formal adjoint} of $D$ is the linear differential operator with constant coefficients given by
\begin{equation*}
  D^\dagger = \sum_{\alpha} (-1)^{\abs{\alpha}} a_{\alpha} \partial^{\alpha}.
\end{equation*}
Here we use the notation $\abs{\alpha} = \alpha_1 + \cdots + \alpha_d$.
The functional $D \eta$ is well defined, because the Schwartz class is stable under the application of
linear differential operators with constant coefficients.
It can be shown that $D \eta$ is again a tempered distribution.
The appearance of $-1$ in $D^{\dagger}$ is justified by the adjoint identity
\begin{equation*}
  \int (Df)(x) g(x) dx = \int f(x) (D^{\dagger}g)(x) dx,
\end{equation*}
which holds for all functions $f,g \in \schwartz{\RR^d}$
and may be obtained via integration by parts.

Naturally, we now define the \emph{distributional derivative} of the field by
\begin{equation*}
  D \varphi(f) = \varphi(D^{\dagger} f) \qquad \forall f \in \schwartz{\RR^4}.
\end{equation*}
As one expects, $D\varphi$ is an operator-valued tempered distribution on $M=\RR^4$.
In terms of creation and annihilation operators we have 
\begin{equation}
  \label{equation:derivative-free-field}
  D \varphi(f) = \frac{1}{\sqrt{2}} \parens*{a(ED^{\dagger}f)^{\dagger} + a(E\overline{D^{\dagger}f})}.
\end{equation}
In Fourier space the operator $D^\dagger$ corresponds to multiplication with the polynomial
\begin{equation}
  \label{equation:d-hat}
  \ft{D}(p) \defequal \sum_{\alpha} i^{\abs{\alpha}} a_{\alpha} (+p^0)^{\alpha_0} (-p^1)^{\alpha_1} (-p^2)^{\alpha_2} (-p^3)^{\alpha_3}
\end{equation}
If $D=\partial^{\mu}$, then $\ft{D}(p) = i @ p_{\!\mu}$, were the potential sign is concealed by lowering the index.

Let $D_1, \ldots, D_r$ be linear differential operators with constant (complex) coefficients.
Plugging~\eqref{equation:derivative-free-field} into the renormalization formula~\eqref{equation:normal-ordering-symmetric} yields
\begin{gather}
  \label{equation:renormalized-product}
  \normord{D_1 \varphi(f) \cdots D_r \varphi(f)} =
  \frac{1}{\sqrt{2^r}}
  \sum_{\sigma \in \SymmetricGroup{r}}
  \sum_{s=0\vphantom{S}}^{r}
  \frac{1}{s!(r-s)!}
  \prod_{i=1\vphantom{S}}^{s} a^\dagger(ED^\dagger_{\sigma(i)}f)
  \prod_{\mathclap{j=s+1\vphantom{S}}}^{r} a(E\overline{D^\dagger_{\sigma(j)}f}).
\end{gather}
%As discussed in \cref{section:quadratic-forms},
%this is a well-defined quadratic form for each fixed test function $f$.

\section{Renormalized Products of the Free Field and~its~Derivatives}

For any given test function $f \in \schwartz{M}$ the renormalized product $\normord{D_1 \varphi(f) \cdots D_r \varphi(f)}$
is well defined as a Fock space operator, but the product is \emph{not} an operator-valued distribution, unless $r=1$.
This is because it has a multi-linear dependence on the test function.
Conceptually, one wishes to treat any physical quantity derived from the quantum field
on the same footing as the field itself.
One construction to obtain an operator-valued distribution,
is described in~\cite{Segal1969}, \cite{Segal1970} and \cite{Klein1973}.
The idea is to take the limit $f \to \delta_x$, where $\delta_x$ is the delta distribution supported at a point $x \in M$,
resulting in the renormalized product at the point $x$, now just a quadratic form,
which is then smeared with \emph{one} instance of $f$.
As we shall see,
this approach incurs significant technical difficulties.

\begin{lemma}{Integral Representation of the Renormalized Product}{renormalized-product-integral-representation}  
  Let $\varphi$ be the free scalar quantum field with mass parameter $m > 0$.
  Let $D_1, \ldots, D_r$ be linear differential operators with constant (complex) coefficients.
  Then, for arbitrary Schwartz functions $f \in \realschwartz{M}$ and Fock states $\psi,\psi' \in \BosonFock{L^2(X_m^+,d\Omega_m)}$,
  we have
  \begin{equation*}
    \innerp{\psi'\!}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi} =
    \int dp_1 \!\cdots dp_r
    \, \ft{f}(p_1) \cdots\! \ft{f}(p_r)
    \, K_{\psi'\!,\psi}(p_1,\ldots,p_r)
  \end{equation*}
  where the \enquote{integral kernel} is given by
  \begin{multline*}
    K_{\psi'\!,\psi}(p_1,\ldots,p_r) =
    \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \sum_{s=0}^{r}
    \delta_{m-s}^{n-(r-s)}
    \ \chi(p_{s+1}) \cdots \chi(p_{r})
    \ P_s(p_1,\ldots,p_r) \\
    \cdot \sqrt{m(m-1) \cdots (m-s+1)} \sqrt{n(n-1) \cdots (n-(r-s)+1)} \\
    \cdot \int dk_1 \cdots dk_{m-s} 
    \ \overline{\psi'_m(k_1,\ldots,k_{m-s},p_1,\ldots,p_s)}
    \ \psi_n(k_1,\ldots,k_{n-(r-s)},p_{s+1},\ldots,p_r)
  \end{multline*}
  \begin{equation*}
    \text{where} \quad \chi(p) = \begin{cases*}
      \overline{\ft{f}(p)} / \ft{f}(p) & if $\ft{f}(p) \ne 0$, \\
      1 & otherwise,
    \end{cases*}
  \end{equation*}
  \begin{multline*}
    \text{and} \quad P_s(p_1,\ldots,p_r) =
    \frac{1}{\sqrt{2^r}}
    \frac{1}{s!(r-s)!}
    \sum_{\sigma \in \SymmetricGroup{r}}
    \ft{D}_{\sigma(1)}(p_1) \cdots \ft{D}_{\sigma(s)}(p_s) \hspace{1.5cm} \\[-1.5ex]
    \cdot \overline{\ft{D}_{\sigma(s+1)}(p_{s+1}) \cdots \ft{D}_{\sigma(r)}(p_r)}.
  \end{multline*}  
\end{lemma}

Note that $K$ has a remaining dependence on $f$ via $\chi$
even though the notation does not indicate this.
This is made explicit in the alternative integral representation
  \begin{equation}
    \label{equation:alternative-integral-representation}
    \begin{multlined}
      \innerp{\psi'\!}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi} =\\
      \hspace{1cm}
      \sum_{s=0}^{r} 
      \int dp_1 \!\cdots dp_r
      \, \ft{f}(p_1) \cdots\! \ft{f}(p_s)
      \, \overline{\ft{f}(p_{s+1}) \cdots\! \ft{f}(p_r)}
      \, \tilde{K}^s_{\psi'\!,\psi}(p_1,\ldots,p_r)
    \end{multlined}
  \end{equation}
  where the integral kernels are given by
  \begin{multline*}
    \tilde{K}^s_{\psi'\!,\psi}(p_1,\ldots,p_r) =
    P_s(p_1,\ldots,p_r)
    \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
    \delta_{m-s}^{n-(r-s)} \\
    \cdot \sqrt{m(m-1) \cdots (m-s+1)} \sqrt{n(n-1) \cdots (n-(r-s)+1)} \\
    \cdot \int dk_1 \cdots dk_{m-s} 
    \ \overline{\psi'_m(k_1,\ldots,k_{m-s},p_1,\ldots,p_s)}
    \ \psi_n(k_1,\ldots,k_{n-(r-s)},p_{s+1},\ldots,p_r).
  \end{multline*}
  This representation will be more convenient for xxx

\begin{myproof}{lemma:renormalized-product-integral-representation}
  From equation~\eqref{equation:renormalized-product},
  applying the definition of the Fock space inner product,
  and moving all creation operators to the left hand side,
  we get
  \begin{multline*}
    \innerp{\psi'}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi} =
    \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
    \frac{1}{\sqrt{2^r}} \sum_{s=0}^{r} \sum_{\sigma \in \SymmetricGroup{r}} \frac{1}{s!(r-s)!} \\
    \cdot \big\langle
    a(ED_{\sigma(1)}^{\dagger}f) \cdots a(ED_{\sigma(s)}^{\dagger}f) \psi_m,
    a(ED_{\sigma(s+1)}^{\dagger}f) \cdots a(ED_{\sigma(r)}^{\dagger}f) \psi_n
    \big\rangle
  \end{multline*}
  Notice that the inner product in the second line
  can only be nonzero if the particle numbers match up
  after the application of the annihilation operators in each argument,
  that is if $m-s=n-(r-s)$.
  By applying~\eqref{equation:multiple-annihilation-operators} twice,
  this expression may be further expanded into
  \begin{gather*}
    \sqrt{m(m-1) \cdots (m-s+1)}
    \sqrt{n(n-1) \cdots (n-(r-s)+1)}
    \int dk_1 \cdots dk_{m-s} \\
    \int dp_1 \cdots dp_s
    \ ED_{\sigma(1)}^{\dagger}f(p_1) \cdots ED_{\sigma(s)}^{\dagger}f(p_s)
    \ \overline{\psi'_m(k_1,\ldots,k_{m-s},p_1,\ldots,p_s)} \\
    \int dp_{s+1} \cdots dp_r
    \ \overline{ED_{\sigma(s+1)}^{\dagger}f(p_{s+1}) \cdots ED_{\sigma(r)}^{\dagger}f(p_r)}
    \ \psi_n(k_1,\ldots,k_{n-(r-s)},p_{s+1},\ldots,p_r)
  \end{gather*}
  Now recall that $E$ stands for Fourier transformation (followed by restriction to the mass shell)
  and that in Fourier space the linear differential operator $D^\dagger$ corresponds to a
  multiplication with the function $\hat{D}$ defined by~\eqref{equation:d-hat}, so that
  \begin{equation*}
    ED_{\sigma(i)}^{\dagger}f(p_i) = \hat{D}_{\sigma(i)}(p_i) \cdot \ft{f} (p_i) 
    \qquad \forall i
  \end{equation*}
  According to Fubini’s Theorem, we may move the $p$-integrals past the $k$-integrals to the very left of the expression.
  This also allows us to move all factors involving $\ft{f}$ or one of the $\ft{D}_{\sigma(i)}$ in front of the $k$-integrals.
  Finally, we move the summation over $s$ to the front.
  By packaging everything not depending on $\ft{f}$ in new functions $\tilde{K}^{s}_{\psi',\psi}$,
  where $s$ runs from $0$ to $r$,
  we arrive at the representation~\eqref{equation:alternative-integral-representation}.

  A single-integral representation is achieved by
  introducing the functions $\chi$s through the substitution $\overline{\ft{f}} = \chi \ft{f}$,
\end{myproof}

In the special case that $D_1 = \cdots = D_r = D$ we have
\begin{equation*}
  P_s(p_1,\ldots,p_r) =
  \frac{1}{\sqrt{2^r}} \parens*{r \atop s\vphantom{y}}
  \ft{D}(p_1) \cdots \ft{D}(p_s)
  \overline{\ft{D}(p_{s+1}) \cdots \ft{D}(p_r)}.
\end{equation*}
In particular, for squares ($r=2$) we have
\begin{equation*}
  P_s(p_1,p_2) = \begin{cases}
    \tfrac{1}{2} \, \overline{\ft{D}(p_1)\ft{D}(p_2)} & s=0 \\
    \phantom{\tfrac{1}{2}} \, \overline{\ft{D}(p_1)}\ft{D}(p_2) & s=1 \\
    \tfrac{1}{2} \, \ft{D}(p_1)\ft{D}(p_2) & s=2
  \end{cases}
\end{equation*}

The following assertion is key to realizing the idea of taking the limit $f \to \delta_x$.

\begin{lemma}{H-bounds for the Integral Kernel}{integral-kernel-h-bound}
  In the setting of \cref{lemma:renormalized-product-integral-representation},
  there exist a constant $C$, and a positive integer $l$,
  such that for arbitrary test functions $f \in \schwartz{M}$
  and states $\psi,\psi' \in \Domain{H^l}$,
  the function $K_{\psi'\!,\psi}$ is integrable (that is, $L^1$)
  and satisfies the $H$-bound
  \begin{equation*}
    \norm{K_{\psi'\!,\psi}}_1 \le
    C \norm{(1+H)^l \psi'} \norm{(1+H)^l \psi}.
  \end{equation*}
  More specifically, it is sufficient to choose $l > rd + r/2$,
  where $d$ is the highest order of differentiation occurring in $D_1, \ldots, D_r$.
\end{lemma}

The Hamilton operator $H$ acts on $n$-particle states $\psi_n$ as follows:
\begin{gather*}
  H \psi_n(p_1,\ldots,p_n) = \parens[\big]{\omega(p_1) + \cdots + \omega(p_n)} \psi_n(p_1,\ldots,p_n) \\
  \shortintertext{where}
  \omega(p) = \omega(\symbfit{p}) = \sqrt{m^2 + \abs{\symbfit{p}}^2} = \sqrt{m^2 + (p^1)^2 + (p^2)^2 + (p^3)^2}.
\end{gather*}
In the following proof it will be convenient to use the abbreviation
\begin{equation*}
  \omega(p_1,\ldots,p_s) \defequal \omega(p_1) + \cdots + \omega(p_n).
\end{equation*}

\begin{myproof}{lemma:integral-kernel-h-bound}
  We have to find an estimate for
  \begin{equation*}
    \norm{K_{\psi'\!,\psi}}_1 =
    \int dp_1 \!\cdots dp_r
    \, \abs{K_{\psi'\!,\psi}(p_1,\ldots,p_r)}.
  \end{equation*}
   We apply the triangle inequalities for sums and integrals
   to the expression for $K_{\psi'\!,\psi}$ given in \cref{lemma:renormalized-product-integral-representation},
   use the fact that $\chi(p)$ has modulus one, make the estimates
   \begin{equation*}
    m(m-1) \cdots (m-s+1) \le m^r
    \quad \text{and} \quad 
    n(n-1) \cdots (n-(r-s)+1) \le n^r,
   \end{equation*}
   and finally reorder the integration with Fubini’s Theorem
   to obtain
  \begin{equation}
    \label{equation:first-estimate}
    \begin{multlined}[c]
      \norm{K_{\psi'\!,\psi}}_1 \le
      \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \sum_{s=0}^{r}
      \delta_{m-s}^{n-(r-s)} \sqrt{m^r n^r} \\
      \hspace{2.5cm} \cdot \int \!dk \int \!dp'\! \int \!dp \,
      \abs*{P_s(p',p) \, \psi'_m(k,p') \, \psi_n(k,p)},
    \end{multlined}
  \end{equation}
  where we have used the abbreviations
  \begin{align*}
    k &= (k_1,\ldots,k_{m-s})
    \quad p' = (p_1,\ldots,p_s)
    \quad p = (p_{s+1},\ldots,p_r) \\
    dk &= dk_1 \cdots dk_{m-s}
    \quad \text{and so on.}
  \end{align*}
   For following discussion we will assume $n-(r-s)=m-s$.

  Observe that $P_s(p_1,\ldots,p_r)$ is a (complex) polynomial
  in the $4r$ variables $p_i^\mu$, $i=1,\ldots,r$, $\mu=0,\ldots,3$.
  Its degree is given by
  \begin{equation*}
    \deg P_s = \sum_{i=1}^r \deg \ft{D}_i
  \end{equation*}
  that is the sum of the highest orders of differentiation
  occurring in each of the operators $D_1, \ldots, D_r$.
  There is no reason to expect arbitrary states $\psi,\psi'$ to temper fast enough
  to counteract this polynomial growth.
  Thus, the integral in~\eqref{equation:first-estimate} will not converge, in general.
  However, if $\psi$ lies in the domain of $H^a$ for some positive integer $a$,
  then we can be sure that $(1+H)^a \psi$ is square integrable, and we have
  \begin{align*}
    \psi_n(k,p) &= \parens[\big]{1+\omega(k,p)} {}^{-a} (1+H)^a \psi_n(k,p) \\
    \psi'_m(k,p') &= \parens[\big]{1+\omega(k,p')} {}^{-a} (1+H)^a \psi'_m(k,p')
  \end{align*}
  We use this to rewrite the integral part of~\eqref{equation:first-estimate} as follows:
  \begin{equation}
    \label{equation:rewritten-integral}
    \int \!dk \int \!dp'\! \int \!dp \, \abs*{F(k,p',p) \, G'(k,p') \, G(k,p)},
  \end{equation}
  where we have introduced the functions
  \begin{align*}
    F(k,p',p) &= \parens[\big]{1+\omega(k,p')} {}^{-a} \parens[\big]{1+\omega(k,p)} {}^{-a} P_s(p',p) \\
    G'(k,p') &= \sqrt{m^r} (1+H)^a \psi'_m(k,p') \\
    G(k,p) &= \sqrt{n^r} (1+H)^a \psi_n(k,p).
  \end{align*}
  Next, we derive an estimate of~\eqref{equation:rewritten-integral}.
  By Cauchy-Schwarz, we have
  \begin{equation*}
    \abs*{\int dp \, \abs{F(k,p',p)G(k,p)}}^2
    \le \int dp \abs{F(k,p',p)}^2
    \cdot \int dp \abs{G(k,p)}^2
  \end{equation*}
  and this implies
  \begin{equation}
    \label{equation:estimate1}
    \int dp' \abs*{\int \!dp \, \abs{F(k,p',p)G(k,p)}}^2
    \le \int dp \abs{G(k,p)}^2
    \cdot \int dp'\! \int dp \abs{F(k,p',p)}^2.
  \end{equation}
  Notice that the second factor is the $L^2$ norm of $F$ with its first argument held fixed:
  \begin{equation}
    \label{equation:norm}
    \norm{F(k,\cdot,\cdot)}_2^2
    = \int dp'\! \int dp \abs{F(k,p',p)}^2
  \end{equation}
  By another application of Cauchy-Schwarz, and using~\eqref{equation:estimate1} and~\eqref{equation:norm}, we obtain
  \begin{align*}
    \text{\eqref{equation:rewritten-integral}} \,
    &= \int \!dk \int \!dp' \abs{G'(k,p')} \int \!dp \, \abs*{F(k,p',p) \, G(k,p)} \\
    &\le \norm{G'}_2 \parens[\bigg]{\int dk \int dp' \abs*{\int \!dp \, \abs*{F(k,p',p)G(k,p)}}^2}^{\frac{1}{2}} \\
    &\le \norm{G'}_2 \parens[\bigg]{\int dk \, \norm{F(k,\cdot,\cdot)}_2^2 \int dp \abs{G(k,p)}^2}^{\frac{1}{2}} \\
    &\le \norm{G'}_2 \norm{G}_2 \sup_{k} \norm{F(k,\cdot,\cdot)}_2
  \end{align*}
  We claim that there exists a positive constant $C_1$ independent of $m$ and $n$ such that
  \begin{equation*}
    \norm{G}_2 \le C_1 \norm{(1+H)^{a+r/2} \psi_n}_2
  \end{equation*}
  and similarly for $G'$.
  Since $\omega$ has a positive lower bound on the mass shell,
  there exists a constant $\epsilon > 0$ independent of $n$ such that
  \begin{equation*}
    1+\omega(k,p) = 1 + \underbrace{\omega(k_1) + \cdots + \omega(k_{m-s}) + \omega(p_{s+1}) + \cdots + \omega(p_r)}_{\text{$(m-s)+(r-s)=n$ terms}} \ge n \epsilon
  \end{equation*}
  Combine this with the fact that
  $N \psi_n = n\psi_n$, where $N$ is the number operator,
  to see
  \begin{equation*}
    \norm{G}_2 = \norm{N^{r/2} (1+H)^a \psi_n}_2 \le \epsilon^{-r/2} \norm{(1+H)^{a+r/2} \psi_n}_2.
  \end{equation*}
  Now set $C_1 = \epsilon^{-r/2}$.
  The proof for $G'$ is analogous.
  In particular, we have shown that both $\norm{G}_{2}$ and $\norm{G'}_{2}$ are finite,
  provided that $\psi_n$ and $\psi'_m$ lie in the domain of $H^l$, where $l=a+r/2$.

  %$H \psi_n(k,p) = (1+\omega(k,p)) \psi_n(k,p)$ 
  %$\norm{(1+H)\psi_n}_2 \ge n \epsilon \norm{\psi_n} = \epsilon \norm{N \psi_n}$ 

  In order to determine conditions for the finiteness of the remaining factor involving $F$,
  it is desirable to have an estimate of the growth of $P_s$ in terms of $\omega(p_1),\ldots,\omega(p_r)$.
  Notice that it is sufficient to make an estimate that is valid on the support of the measure $\Omega_m$, that is, the mass shell $X_m^+$,
  since $F$ appears in an integral with respect to $p_1,\ldots,p_r$.
  For an arbitrary point $q$ on the mass shell $X_m^+$ we have
  \begin{align*}
    q^{0} &= \omega(q) \\
    \abs{q^{\mu}} &\le \omega(q) \quad \mu = 1,2,3.
  \end{align*}
  Moreover, $\omega(q)$ has a positive lower bound on $X_m^+$, so that
  for all exponents $a,b \in \NN$ with $a < b$ there
  exists a constant $c_{a,b}$ such that $\omega(q)^a \le c_{a,b}\, \omega(q)^b$.
  This allows us to make the estimate
  \begin{equation}
    \label{equation:polynomial-estimate}
    \abs{P_s(p_1,\ldots,p_r)} \le C_s \prod_{i=1}^r \omega(p_i)^{d_i} \quad \text{where}\ d_i = \deg \ft{D}_i,
  \end{equation}
  and $C_s$ is a constant independent of $m$ and $n$.
  By the Arithmetic Mean-Geometric Mean Inequality, we have
  \begin{gather}
    \sqrt[s]{\omega(p_1) \cdots \omega(p_s)} 
    \le \frac{\omega(p_1) + \cdots + \omega(p_s)}{s} 
    \le \omega(p') \le 1 + \omega(k,p'), \nonumber\\
    \shortintertext{hence}
    \label{equation:one-plus-omega-estimate1}
    \parens[\big]{1+\omega(k,p')} {}^{-a}
    \le \parens[\big]{\omega(p_1) \cdots \omega(p_s)} {}^{-a/s}, \\
    \shortintertext{and similarly}
    \label{equation:one-plus-omega-estimate2}
    \parens[\big]{1+\omega(k,p)} {}^{-a}
    \le \parens[\big]{\omega(p_1) \cdots \omega(p_{r-s})} {}^{-a/(r-s)}.
  \end{gather}
  The estimates~\eqref{equation:polynomial-estimate},~\eqref{equation:one-plus-omega-estimate1} and~\eqref{equation:one-plus-omega-estimate2} entail
  \begin{equation*}
    \abs{F(k,p',p)} \le C_s
    \prod_{i=1}^{s} \omega(p_i)^{d_i-a/s} 
    \prod_{j=s+1}^{r} \omega(p_j)^{d_j-a/(r-s)}.
  \end{equation*}
  Since the right hand side does not depend on $k$, and the $p$-variables are separated,
  the problem reduces to proving that
  \begin{equation}
    \label{equation:integral-finite}
    \int \omega(q)^{-2b} \,d\Omega(q) < \infty
  \end{equation}
  for $b$ large enough.
  Recall that $d \Omega(q) = \omega(q)^{-1} d^3 \symbfit{q}$.
  By transformation to spherical coordinates, we find that~\eqref{equation:integral-finite}
  is equivalent to
  \begin{equation*}
    \int \frac{r^2}{(m^2 + r^2)^{b+1/2}} \,dr < \infty
  \end{equation*}
  It it well known that this holds for $b > 1$.
  $a > r d$ 

  \begin{equation}
    \label{equation:intermediate-result}
    \norm{K_{\psi'\!,\psi}}_1 \le
    \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
    \sum_{s=0}^{r} \delta_{m-s}^{n-(r-s)} C_s
    \underbracket{\norm{(1+H)^l \psi'_m}_2}_{a'_m}
    \underbracket{\norm{(1+H)^l \psi_n}_2}_{a_n}
  \end{equation}
  We introduce auxiliary variables $a'_m, a_n$ as shown above, and for convenience set $a_n = 0$ whenever $n < 0$.
  Using this, we rewrite the right hand side of~\eqref{equation:intermediate-result}
  and apply the Cauchy-Schwarz Inequality for sequences as follows:
  \begin{equation*}
    \sum_{s=0}^{r} C_s \sum_{m = 0}^{\infty} a'_m \cdot a_{m+r-2s} \le
    \sum_{s=0}^{r} C_s \sqrt{\sum_{m=0}^{\infty} a'^2_m \sum_{n=0}^{\infty} a^2_n}
  \end{equation*}
  To complete the proof, observe that
  \begin{equation*}
    \norm{(1+H)^l \psi'}_2 =
    \sqrt{\sum_{m=0}^{\infty} \norm{(1+H)^l \psi'_m}_2^2} =
    \sqrt{\sum_{m=0}^{\infty} a'^2_m},
  \end{equation*}
  and similar for $\psi$, by definition of the inner product
  and because $((1+H)^l \psi')_m = (1+H)^l \psi'_m$ for all $m$.
\end{myproof}

Now we are in a position to implement the idea of taking the limit $f \to \delta_x$.
Recall that any tempered distribution, and in particular Dirac distributions, may be approximated by Schwarz test functions.

\begin{lemma}{Renormalized Product at a Point}{}
  In the setting of \cref{lemma:renormalized-product-integral-representation},
  assume that $\psi,\psi'$ are in $\Domain{H^l}$.
  Let $x$ be any point in $M$ and let $\delta_x \in \tempdistrib{M}$ be the Dirac distribution supported at $x$.
  Then the limit
  \begin{equation*}
    \lim_{f \to \delta_x}
    \innerp{\psi'\!}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi}
  \end{equation*}
  exists and depends continuously on $x$.
\end{lemma}

\begin{proof}
  Since the Fourier transformation of tempered distribution
  is a continuous mapping $\tempdistribnoarg \to \tempdistribnoarg$,
  we have $\ft{f} \to \FT{\delta_x}$ whenever $f \to \delta_x$ in the topology of $\tempdistribnoarg$.
  Consequently, $\abs{\ft{f}}$ remains bounded by some constant $C$ while taking the limit.

  According to \cref{lemma:renormalized-product-integral-representation} we have
  \begin{equation*}
    \innerp{\psi'\!}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi}
    = \int dp_1 \!\cdots dp_r
    \, \ft{f}(p_1) \cdots\! \ft{f}(p_r)
    \, K_{\psi'\!,\psi}(p_1,\ldots,p_r)
  \end{equation*}
  The integrand is dominated by the function $C^r\abs{K_{\psi'\!,\psi}(p_1,\ldots,p_r)}$,
  which has finite integral by \cref{lemma:integral-kernel-h-bound}.

  %Moreover, the integrand converges pointwise to $K_{\psi'\!,\psi}(p_1,\ldots,p_r)$, since $\ft{f} \to 1$ when $f \to \delta_x$.
  %\todo{With of choice of FT constants, $\ft{f} \to 1/(2\pi)^2$. Change here or change def?}

  Recall that $\ft{\delta} = 1$, and thus $\FT{\delta_x}(p) = e^{ix \cdot p}$ for all $p \in M$.
  This shows that the integrand converges pointwise to the function
  \begin{equation*}
    F(p_1,\ldots,p_r) =
    \sum_{s=0}^r
    e^{ix \cdot (p_1 + \cdots + p_s)}
    e^{-ix \cdot (p_{s+1} + \cdots + p_r)}
    \tilde{K}^s_{\psi'\!,\psi}(p_1,\ldots,p_r)
  \end{equation*}
  The Dominated Convergence Theorem implies that $F$ is integrable and that
  \begin{equation}
    \label{equation:intx}
    \lim_{f \to \delta_x}
    \innerp{\psi'\!}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi}
    = \int dp_1 \!\cdots dp_r \, F(p_1,\ldots,p_r).
  \end{equation}
  In particular, the limit exists.
\end{proof}

\begin{definition}{Renormalized Product as a QF-valued distribution}{}
  Adopt the assumptions of the foregoing Lemma.
  We define two mappings that intentionally share the same name by
  \begin{gather*}
    \normord{D_1 \varphi \cdots D_r \varphi} \ \vcentcolon \ 
    M \to \QF{\mathcal{D}(H^l)} \\
    \innerp{\psi'\!}{\normord{D_1 \varphi \cdots D_r \varphi}(x) \,\psi}
    = \lim_{f \to \delta_x}
    \innerp{\psi'\!}{\normord{D_1 \varphi(f) \cdots D_r \varphi(f)} \,\psi}
  \end{gather*}
  and
  \begin{equation*}
    \normord{D_1 \varphi \cdots D_r \varphi} \ \vcentcolon \ 
    \schwartz{M} \to \QF{\mathcal{D}(H^l)}
  \end{equation*}
  \begin{equation*}
    \innerp{\psi'\!}{\normord{D_1 \varphi \cdots D_r \varphi}(f) \,\psi} =
    \int_M \!dx \ f(x) \ \innerp{\psi'\!}{\normord{D_1 \varphi \cdots D_r \varphi}(x) \,\psi}
  \end{equation*}
\end{definition}

\begin{proposition}{}{}  
  Let $\varphi$ be a free quantum field.
  Let $D_1, \ldots, D_r$ be linear differential operators with constant (complex) coefficients.
  Suppose that $l$ is a large enough positive integer.
  Then we have for all states $\psi,\psi' \in \Domain{H^l}$
  \begin{multline*}
    \innerp{\psi'\!}{\normord{D_1 \varphi \cdots D_r \varphi}(f) \,\psi} = \\
    = \int dp_1 \!\cdots dp_r
    \sum_{s=0}^{r}
    \, \ft{f}(p_1 + \cdots + p_s - p_{s+1} - \cdots - p_r)
    \, \tilde{K}_{\psi'\!,\psi}^{s}(p_1,\ldots,p_r)
  \end{multline*}
  where
  \begin{multline*}
    \tilde{K}_{\psi'\!,\psi}^{s}(p_1,\ldots,p_r) =
    P_s(p_1,\ldots,p_r)
    \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
    \delta_{m-s}^{n-(r-s)} \\
    \cdot \sqrt{m(m-1) \cdots (m-s+1)} \sqrt{n(n-1) \cdots (n-(r-s)+1)} \\
    \cdot \int dk_1 \cdots dk_{m-s} 
    \ \overline{\psi'_m(k_1,\ldots,k_{m-s},p_1,\ldots,p_s)}
    \ \psi_n(k_1,\ldots,k_{n-(r-s)},p_{s+1},\ldots,p_r)
  \end{multline*}
  and $P_s(p_1,\ldots,p_r)$ is defined as before.
\end{proposition}

\begin{proof}
  This follows directly from the definition and~\eqref{equation:intx}.
\end{proof}


\section{Definition of the Stress Tensor}

In the theory of a real scalar field $\phi$ of mass $m$,
the Lagrangian density of the Klein-Gordon action is given by
\begin{equation}
  \label{lagrangian-density}
  \mathcal{L} = \frac{1}{2} \parens{\partial^{\mu} \phi \partial_{\mu} \phi - m^2 \phi^2} 
\end{equation}
and the \emph{canonical stress-energy tensor} is defined by
\begin{equation*}
  T^{\mu}_{\nu} = \frac{\partial \mathcal{L}}{\partial(\partial_{\mu} \phi)} \partial_{\nu} \phi - \delta^{\mu}_{\nu} \mathcal{L}
\end{equation*}
Raising the index $\nu$ and inserting \cref{lagrangian-density} yields
\begin{equation*}
  T^{\mu\nu} = \partial^{\mu}\phi \partial^{\nu}\phi + \frac{1}{2} \eta^{\mu\nu} \parens*{m^2 \phi^2 - \partial_{\lambda}\phi \partial^{\lambda}\phi} 
\end{equation*}
The \emph{energy density}:
\begin{equation*}
  \energydensity = T^{00} = \frac{1}{2} \parens*{\sum_{\mu=0}^{3} (\partial^{\mu}\phi)^2  + m^2 \phi^2} 
\end{equation*}
The discussion in the previous section enables us to define
the \emph{renormalized stress-energy tensor} of a free scalar field $\varphi$ by
\begin{equation*}
  T^{\mu\nu} = \normord{\partial^{\mu}\varphi \partial^{\nu}\varphi + \frac{1}{2} \eta^{\mu\nu} \parens*{m^2 \varphi^2 - \partial_{\lambda}\varphi \partial^{\lambda}\varphi}}
\end{equation*}
and this is a quadratic form.
In particular, the energy density is
\begin{equation*}
  \energydensity  = \frac{1}{2} \sum_{\mu=0}^{3} \normord{(\partial^{\mu}\varphi)^2}  + \frac{1}{2} m^2 \normord{\varphi^2}
\end{equation*}

Now let $f$ be a real-valued Schwarz function,
and let $\psi,\psi' \in \Domain{H^l}$ for a large enough integer $l$.
For convenience we introduce the functions
\begin{multline*}
  L_{\psi'\!,\psi}^{s}(p_1,\ldots,p_r) =
  \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
  \delta_{m-s}^{n-(r-s)} \\
  \cdot \sqrt{m(m-1) \cdots (m-s+1)} \sqrt{n(n-1) \cdots (n-(r-s)+1)} \\
  \cdot \int dk_1 \cdots dk_{m-s} 
  \ \overline{\psi'_m(k_1,\ldots,k_{m-s},p_1,\ldots,p_s)}
  \ \psi_n(k_1,\ldots,k_{n-(r-s)},p_{s+1},\ldots,p_r)
\end{multline*}
so that
\begin{equation*}
  \tilde{K}_{\psi'\!,\psi}^{s}(p_1,\ldots,p_r) =
  P_s(p_1,\ldots,p_r) \,
  L_{\psi'\!,\psi}^{s}(p_1,\ldots,p_r).
\end{equation*}
Since the energy density only contains squares, it suffices to consider $r=2$.

\begin{multline*}
  \innerp{\psi'\!}{\normord{\varphi^2}(f) \,\psi} = \int dp_1 dp_2 \,
  \ft{f}(-p_1 - p_2) \, \tfrac{1}{2} L_{\psi'\!,\psi}^{0}(p_1,p_2) \\
  + \ft{f}(p_1 - p_2) \, L_{\psi'\!,\psi}^{1}(p_1,p_2)
  + \ft{f}(p_1 + p_2) \, \tfrac{1}{2} L_{\psi'\!,\psi}^{2}(p_1,p_2) \\
  \innerp{\psi'\!}{\normord{(\partial^{\mu}\varphi)^2}(f) \,\psi} = \int dp_1 dp_2 \,
  \ft{f}(-p_1 - p_2) \, \parens[\big]{-\tfrac{1}{2} (p_1)_{\mu} (p_2)_{\mu}} \, L_{\psi'\!,\psi}^{0}(p_1,p_2) \\
  + \ft{f}(p_1 - p_2) \, \parens[\big]{+(p_1)_{\mu} (p_2)_{\mu}} \, L_{\psi'\!,\psi}^{1}(p_1,p_2)
  + \ft{f}(p_1 + p_2) \, \parens[\big]{-\tfrac{1}{2} (p_1)_{\mu} (p_2)_{\mu}} \, L_{\psi'\!,\psi}^{2}(p_1,p_2)
\end{multline*}

\begin{equation*}
  \bar{p} \defequal \eta p = (p^0,-\Vector{p}) 
\end{equation*}

\begin{equation*}
  \sum_{\mu = 0}^{3} p_{\mu} p'_{\mu} = \bar{p} \cdot p' 
\end{equation*}

\begin{proposition}{}{}
  \begin{multline*}
    \innerp{\psi'\!}{\energydensity(f) \, \psi} =
    \frac{1}{4} \int dp \, dp'
    (m^2 + \bar{p} \cdot p')
    \bracks[\big]{2 \ft{f}(p - p') L^1_{\psi'\!,\psi}(p,p')} \\
    + (m^2 -\bar{p} \cdot p')
    \bracks[\big]{\ft{f}(- p - p') L^0_{\psi'\!,\psi}(p,p') + \ft{f}(p + p') L^2_{\psi'\!,\psi}(p,p')} 
  \end{multline*}
\end{proposition}

\begin{theorem}{TODO}{}  
  Let $\varphi$ be a free quantum field.
  Let $D_1, \ldots, D_r$ be linear differential operators with constant (complex) coefficients.
  Then we have for all test functions $f \in \schwartz{M}$
  \begin{multline*}
    \normord{D_1 \varphi \cdots D_r \varphi}(f) \QFequal \int dp_1 \!\cdots dp_r \sum_{s=0}^{r}
    P_s(p_1,\ldots,p_r) \, \ft{f}(p_1 + \cdots + p_s - p_{s+1} - \cdots - p_r) \\
    \cdot a^\dagger(p_1) \cdots a^\dagger(p_s) a(p_{s+1}) \cdots a(p_r) 
  \end{multline*}
  as quadratic forms on $\Domain{H^l}$, where
  \begin{multline*}
    \quad P_s(p_1,\ldots,p_r) =
    \frac{1}{\sqrt{2^r}}
    \frac{1}{s!(r-s)!}
    \sum_{\sigma \in \SymmetricGroup{r}}
    \ft{D}_{\sigma(1)}(p_1) \cdots \ft{D}_{\sigma(s)}(p_s) \hspace{1.5cm} \\[-1.5ex]
    \cdot \overline{\ft{D}_{\sigma(s+1)}(p_{s+1}) \cdots \ft{D}_{\sigma(r)}(p_r)}.
  \end{multline*}  
\end{theorem}

\begin{proposition}{}{energy-density}
  \begin{multline*}
    \energydensity(f) \QFequal \frac{1}{4} \int dp \, dp' \,
    (m^2 + \bar{p} \cdot p')
    \bracks[\Big]{2\ft{f}(p-p') a^\dagger(p) a(p')} + {} \\
    + (m^2 - \bar{p} \cdot p')
    \bracks[\Big]{\ft{f}(p+p') a(p) a(p') + \ft{f}(-p-p') a^\dagger(p) a^\dagger(p')} 
  \end{multline*}
\end{proposition}

Observe that $\bar{p} \cdot p' = p^0p'^0 + p^1p'^1 + p^2p'^2 + p^3p'^3$ is symmetric in $p$ and $p'$.
Consequently, we could rewrite the first bracketed expression as $2\ft{f}(p-p') a^\dagger(p) a(p')$ 

\begin{proposition}{}{}
  The Fock vacuum $\FockVacuum$ lies in the domain of $\energydensity(f)\QFop{}$
  for all test functions $f \in \schwartz{M}$
  and $\energydensity(f)\QFop{}\FockVacuum$ is the vector $\psi$ defined by
  \begin{equation*}
    \psi_2(p,p') = \frac{\sqrt{2}}{4} (m^2 - \bar{p} \cdot p') \ft{f}(-p-p')
  \end{equation*}
  and $\psi_n \equiv 0$ for $n \ne 2$.
\end{proposition}

\begin{proof}
  Since $\ft{f}$ is a Schwarz function, there exists for each $N \in \NN$  a positive constant $C_N$ such that
  \begin{equation*}
    \abs{\ft{f}(q)}^2 \le \frac{C_N}{1+(q^0)^{2N}}
  \end{equation*}
  For $p,p'$ on the mass shell we have
  \begin{equation*}
    m^2 - \bar{p} \cdot p' =
    m^2 - \sqrt{m^2 + \norm{\Vector{p}}^2} \sqrt{m^2 + \norm{\Vector{p'}}^2} + \Vector{p} \cdot \Vector{p}',
  \end{equation*}
  where on the right hand side \enquote*{$\cdot$} stands for the Euclidean scalar product for three-vectors.
  The Cauchy-Schwarz Inequality (applied twice) shows
  \begin{equation*}
    m^2 + \Vector{p} \cdot \Vector{p}'
    \le m^2 + \norm{\Vector{p}} \norm{\Vector{p}'}
    \le \sqrt{m^2 + \norm{\Vector{p}}^2} \sqrt{m^2 + \norm{\Vector{p'}}^2}.
  \end{equation*}
  \begin{equation*}
    \abs{m^2 - \bar{p} \cdot p'}
    \le 2 \sqrt{m^2 + \norm{\Vector{p}}^2} \sqrt{m^2 + \norm{\Vector{p'}}^2}
  \end{equation*}
\end{proof}

\section{Essential Selfadjointness of Renormalized Products}

\begin{lemma}{H-Bounds for the Renormalized Product}{}
  \begin{equation*}
    \abs{\innerp{\psi'\!}{\normord{D_1 \varphi \cdots D_r \varphi}(f)\psi}} \le
    C \norm{(I+H)^l \psi} \norm{(I+H)^l \psi'}
  \end{equation*}
\end{lemma}

\begin{proof}
  This proof is nearly identical to that of \cref{lemma:integral-kernel-h-bound}
  and we will only cover the differences.
  \begin{equation*}
    \begin{multlined}[c]
    \abs{\innerp{\psi'\!}{\normord{D_1 \varphi \cdots D_r \varphi}(f)\psi}} \le
      \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \sum_{s=0}^{r}
      \delta_{m-s}^{n-(r-s)} \\
      \hspace{2.5cm} \cdot
    \int \!dk \int \!dp'\! \int \!dp \, \abs*{F(k,p',p) \, G'(k,p') \, G(k,p)},
    \end{multlined}
  \end{equation*}
  where
  \begin{multline}
    F(k,p',p) = \parens[\big]{1+\omega(k,p')} {}^{-a} \parens[\big]{1+\omega(k,p)} {}^{-a} P_s(p',p) \\
    \cdot \ft{f}(p_1 + \cdots + p_s - p_{s+1} - \cdots - p_r)
  \end{multline}
  and $G$ and $G'$ are defined as before.
  All we have to do, is verifying that
  \begin{equation*}
    \sup_k \norm{F(k,\cdot,\cdot)}_2 < \infty
  \end{equation*}
  for a sufficiently large integer $a$.
  Then we obtain the desired $H$-bound with $l=a+r/2$.

  Recall that the Schwartz class is preserved by Fourier transform, translation and multiplication with polynomials.
  Moreover, it is well known that Schwartz functions are square-integrable with respect to the Lorentz invariant measure on the mass shell.
  Hence,
  \begin{equation*}
    \int dp_1 \abs{\ft{f}(p_1 + \cdots + p_s - p_{s+1} - \cdots - p_r)}^2
  \end{equation*}
  is bounded by a constant independent of $p_2,\ldots,p_r$.
\end{proof}

\section{Covariance}

\begin{equation*}
  f_g(x) \defequal f(g^{-1} x) \qquad
  x \in M, \quad g \in \ProperOrthochronousPoincareGroup.
\end{equation*}

\begin{theorem}{Covariance}{covariance-renormalized-product}
  \begin{equation*}
    U(g) \,\normord{D_1 \varphi \cdots D_r \varphi}(f)\, U(g)^{-1}
    = \normord{D_1 \varphi \cdots D_r \varphi}(f_g)
  \end{equation*}
\end{theorem}

% vim: syntax=mytex