# Optimal density compensation factors for the reconstruction of the Fourier
transform of bandlimited functions

Melanie Kircheis111Corresponding author<EMAIL_ADDRESS>Chemnitz University of Technology, Faculty of Mathematics, D–09107 Chemnitz,
Germany Daniel<EMAIL_ADDRESS>Chemnitz University of
Technology, Faculty of Mathematics, D–09107 Chemnitz, Germany

###### Abstract

An inverse nonequispaced fast Fourier transform (iNFFT) is a fast algorithm to
compute the Fourier coefficients of a trigonometric polynomial from
nonequispaced sampling data. However, various applications such as magnetic
resonance imaging (MRI) are concerned with the analogous problem for
bandlimited functions, i. e., the reconstruction of point evaluations of the
Fourier transform from given measurements of the bandlimited function. In this
paper, we review an approach yielding exact reconstruction for trigonometric
polynomials up to a certain degree, and extend this technique to the setting
of bandlimited functions. Here we especially focus on methods computing a
diagonal matrix of weights needed for sampling density compensation.

_Key words_ : inverse nonequispaced fast Fourier transform, nonuniform fast
Fourier transform, direct inversion, density compensation, bandlimited
functions, iNFFT, NFFT, NUFFT

AMS _Subject Classifications_ :  65Txx, 65T50, 65F05.

##  1 Introduction

The nonequispaced fast Fourier transform (NFFT) is a fast algorithm to
evaluate a trigonometric polynomial

$\displaystyle
f(\boldsymbol{x})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
(1.1)

with given Fourier coefficients $\hat{f}_{\boldsymbol{k}}\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, at nonequispaced points
$\boldsymbol{x}_{j}$, $j=1,\dots,N$, with $N\in\mathbb{N}$, where for
$\boldsymbol{M}\coloneqq(M,\dots,M)^{T}$, $M\in 2\mathbb{N}$, we define the
multi-index set
$\mathcal{I}_{\boldsymbol{M}}\coloneqq\mathbb{Z}^{d}\cap\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}$
with cardinality $|\mathcal{I}_{\boldsymbol{M}}|=M^{d}$. For more information
see [8, pp. 377-381] and references therein. In case of equispaced points
$\boldsymbol{x}_{j}$ and $|\mathcal{I}_{\boldsymbol{M}}|=N$, this evaluation
can be realized by means of the well-known fast Fourier transform (FFT); an
algorithm that is invertible. However, various applications need to perform an
inverse nonequispaced fast Fourier transform (iNFFT), i. e., compute the
Fourier coefficients $\hat{f}_{\boldsymbol{k}}$ from given function
evaluations $f(\boldsymbol{x}_{j})$ of the trigonometric polynomial (1.1).
Hence, we are interested in an inversion also for nonequispaced data. In
contrast to iterated procedures, where multiple applications of the NFFT are
needed to compute a solution, we consider so-called direct methods, where for
a fixed set of points $\boldsymbol{x}_{j}$, $j=1,\dots,N$, the reconstruction
can be realized with the same number of arithmetic operations as a single
application of an adjoint NFFT. To achieve this, a certain precomputational
step is compulsory, since the adjoint NFFT does not yield an inversion of the
NFFT by itself, see [4, Sec. 3]. Although this precomputations might be rather
costly, they need to be done only once for a given set of points
$\boldsymbol{x}_{j}$, $j=1,\dots,N$, while the actual reconstruction step is
very efficient. Thus, direct methods are especially beneficial in case of
multiple applications with fixed points.

For the setting of trigonometric polynomials several approaches are known,
while in this paper we focus on the analogous problem for bandlimited
functions with

$\displaystyle
f(\boldsymbol{x})=\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v},\quad\boldsymbol{x}\in\mathbb{R}^{d}.$
(1.2)

In various applications such as magnetic resonance imaging (MRI), cf. [1], the
aim is the reconstruction of point evaluations
$\hat{f}(\boldsymbol{k})\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, of an object $\hat{f}$ from
given measurements $f(\boldsymbol{x}_{j})$, $j=1,\dots,N$, of the form (1.2).
To this end, known approaches for the direct inversion of the NFFT shall be
extended to this setting. In particular, we investigate methods computing a
diagonal matrix of weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N$, which are
needed for sampling density compensation due to the nonequispaced sampling
$\boldsymbol{x}_{j}$, $j=1,\dots,N$.

Therefore, this paper is organized as follows. Firstly, in Section 2 we review
the approach in [4, Sec. 3.1], that leads to an exact reconstruction for all
trigonometric polynomials (1.1) of degree $\boldsymbol{M}$. Subsequently, in
Section 3 we aim for an analogous approach for bandlimited functions $f\in
L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ as well. Finally, in Section
4 we show some numerical examples to investigate the accuracy of our
approaches.

##  2 Trigonometric polynomials

Let $\mathbb{T}^{d}\coloneqq\mathbb{R}^{d}\setminus\mathbb{Z}^{d}$ denote the
$d$-dimensional torus with $d\in\mathbb{N}$. The inner product of two vectors
shall be defined as usual as $\boldsymbol{k}\boldsymbol{x}\coloneqq
k_{1}x_{1}+\dots+k_{d}x_{d}$. For the reconstruction of the Fourier
coefficients $\hat{f}_{\boldsymbol{k}}\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, from function evaluations
$f(\boldsymbol{x}_{j})$, $j=1,\dots,N$, of trigonometric polynomials $f\in
L_{2}(\mathbb{T}^{d})$ of degree $\boldsymbol{M}$, see (1.1), the following
result is known, cf. [4, Cor. 3.4].

######  Theorem 2.1.

Let $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$ and
$\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N$, be given. Then the
density compensation factors $w_{j}\in\mathbb{C}$ satisfying

$\displaystyle\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\delta_{\boldsymbol{0},\boldsymbol{k}},\quad\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{2M}}},$
(2.1)

with the Kronecker delta $\delta_{\boldsymbol{\ell},\boldsymbol{k}}$, are
optimal, since for all trigonometric polynomials of degree $\boldsymbol{M}$,
see (1.1), an exact reconstruction of the Fourier coefficients
$\hat{f}_{\boldsymbol{k}}$ is given by

$\displaystyle\hat{f}_{\boldsymbol{k}}=h_{\boldsymbol{k}}^{\mathrm{w}}\coloneqq\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$
(2.2)

###### Proof.

It is known that
$\\{\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{\ell}\boldsymbol{x}}\colon\boldsymbol{\ell}\in\mathbb{Z}^{d}\\}$
forms an orthonormal basis of the Hilbert space $L_{2}(\mathbb{T}^{d})$ of all
$1$-periodic, complex-valued functions, see [8, p. 161]. To achieve an exact
reconstruction for all trigonometric polynomials (1.1) with maximum degree
$\boldsymbol{M}$ it suffices to consider the set of basis functions with
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$. For each of these basis
functions with fixed $\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$ we
have

$\displaystyle
h_{\boldsymbol{k}}^{\mathrm{w}}=\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$

and

$\displaystyle\hat{f}_{\boldsymbol{k}}$
$\displaystyle=c_{\boldsymbol{k}}(f)=\int\limits_{\mathbb{T}^{d}}\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}}\,\mathrm{d}{\boldsymbol{x}}=\delta_{\boldsymbol{\ell},\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$

where $\delta_{\boldsymbol{\ell},\boldsymbol{k}}$ denotes the Kronecker delta.
Thus, to obtain (2.2) we need to assure that the weights satisfy

$\displaystyle\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}=\delta_{\boldsymbol{\ell},\boldsymbol{k}},\quad\boldsymbol{\ell},\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$
(2.3)

Since for $\boldsymbol{k},\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$ we
have $(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{2M}}$, the
property (2.3) is fulfilled by (2.1). In other words, by defining the
nonequispaced Fourier matrix

$\displaystyle\boldsymbol{A}=\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{M}}|}\coloneqq\left(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\right)_{j=1,\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}^{N}\
\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ (2.4)

as well as $\boldsymbol{w}\coloneqq(w_{j})_{j=1}^{N}$, an exact solution to
the linear system of equations

$\displaystyle\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\boldsymbol{e}_{\boldsymbol{0}}\coloneqq\left(\delta_{\boldsymbol{0},\boldsymbol{k}}\right)_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{2M}}}}$
(2.5)

yields an exact reconstruction for all trigonometric polynomials (1.1) of
degree $\boldsymbol{M}$. As already mentioned in [3, Sec. 3.1] an exact
solution to (2.5) can only be found if (2.5) is an underdetermined system with
$|\mathcal{I}_{\boldsymbol{2M}}|\leq N$. ∎

Thus, in case $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$ and the condition number
of $\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}$ is sufficiently small,
the optimal density compensation factors $\boldsymbol{w}\in\mathbb{C}^{N}$ can
be found as solution to the linear system of equations (2.5). As mentioned in
[4, Sec. 3.1.1] the unique solution is then given by the normal equations of
second kind

$\displaystyle\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{v}=\boldsymbol{e}_{\boldsymbol{0}},\quad\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{v}=\boldsymbol{w},$
(2.6)

and can efficiently be computed by an iteration procedure combining the CG
algorithm and the NFFT.

However, when $|\mathcal{I}_{\boldsymbol{2M}}|>N$ there is no theoretical
guaranty. In addition, the numerics in [4] have shown that the least squares
solution satisfying the normal equations of first kind

$\displaystyle\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{e}_{\boldsymbol{0}},$
(2.7)

is not a good approximation. This is why, we recommend another computation
scheme in this setting. To this end, note that (2.5) implies by (2.3) that

$\displaystyle\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
(2.8)

with the weight matrix
$\boldsymbol{W}\coloneqq\mathrm{diag}(w_{j})_{j=1}^{N}$. Hence, in case a
solution to (2.5) cannot be found using (2.6), one may try to find suitable
weights utilizing (2.8) by minimizing the Frobenius norm
$\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}$.
As shown in [4, Sec. 3.4.2] this minimizer can be obtained by solving
$\boldsymbol{S}\boldsymbol{w}=\boldsymbol{b}$, where

$\displaystyle\boldsymbol{S}$
$\displaystyle\coloneqq\Big{(}\left|\left[\boldsymbol{A}\boldsymbol{A}^{*}\right]_{j,s}\right|^{2}\Big{)}_{j,s=1}^{N}\quad\text{
and
}\quad\boldsymbol{b}=|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{1}_{N}.$
(2.9)

Note that by means of the definition (2.4) as well as
$\boldsymbol{f}\coloneqq\left(f(\boldsymbol{x}_{j})\right)_{j=1}^{N}$ and
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
the reconstruction (2.2) can be denoted as
$\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$. For
further details see [4].

##  3 Bandlimited functions

Since various applications are concerned with bandlimited functions instead of
trigonometric polynomials, we now aim to extend the method from Section 2 to
bandlimited functions $f\in L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$
with maximum bandwidth $\boldsymbol{M}$, i. e., functions whose (continuous)
Fourier transform

$\displaystyle\hat{f}(\boldsymbol{v})\coloneqq\int\limits_{\mathbb{R}^{d}}f(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x},\quad\boldsymbol{v}\in\mathbb{R}^{d},$
(3.1)

is supported on $\left[-\frac{M}{2},\frac{M}{2}\right)^{d}$. Thus, their
inverse Fourier transform

$\displaystyle
f(\boldsymbol{x})=\int\limits_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v},\quad\boldsymbol{x}\in\mathbb{R}^{d},$

can be written as in (1.2).

###  3.1 Reconsideration as trigonometric polynomials

To find a suitable reconstruction technique for a bandlimited function $f$, we
consider its 1-periodized version

$\displaystyle\tilde{f}(\boldsymbol{x})\coloneqq\sum_{\boldsymbol{r}\in\mathbb{Z}^{d}}f(\boldsymbol{x}+\boldsymbol{r}).$

Note that $\tilde{f}\in L_{2}(\mathbb{T}^{d})$ is uniquely representable in
form of its absolute convergent Fourier series

$\displaystyle\tilde{f}(\boldsymbol{x})\coloneqq\sum_{\boldsymbol{k}\in\mathbb{Z}^{d}}c_{\boldsymbol{k}}(\tilde{f})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}},$
(3.2)

see [8, Thm. 4.5], where the Fourier coefficients are given by

$\displaystyle c_{\boldsymbol{k}}(\tilde{f})$
$\displaystyle=\int\limits_{\mathbb{T}^{d}}\tilde{f}(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}=\int\limits_{\mathbb{R}^{d}}f(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}=\hat{f}(\boldsymbol{k}),\quad\boldsymbol{k}\in\mathbb{Z}^{d},$
(3.3)

cf. (3.1). Moreover, it is known that $f$ is a bandlimited function with
bandwidth $\boldsymbol{M}$, i. e., we have $\hat{f}(\boldsymbol{k})=0$,
$\boldsymbol{k}\in\mathbb{Z}^{d}\setminus\mathcal{I}_{\boldsymbol{M}}$.
Therefore, the periodic function (3.2) in fact is a trigonometric polynomial
of degree $\boldsymbol{M}$ as in (1.1), which makes it reasonable to utilize
the result from Section 2 to reconstruct the Fourier coefficients (3.3).
Namely, an exact solution to the linear system (2.5) yields an exact
reconstruction

$\displaystyle\hat{f}(\boldsymbol{k})=c_{\boldsymbol{k}}(\tilde{f})=\sum_{j=1}^{N}w_{j}\,\tilde{f}(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$
(3.4)

Note that using the nonequispaced Fourier matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in
(2.4), the weight matrix $\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}$ as
well as the vectors
$\boldsymbol{\tilde{f}}\coloneqq(\tilde{f}(\boldsymbol{x}_{j}))_{j=1}^{N}$ and
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
the reconstruction (3.4) can be denoted as
$\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{\tilde{f}}$.
In addition, by (3.2) we have
$\boldsymbol{\tilde{f}}=\boldsymbol{A}\boldsymbol{\hat{f}}$, such that (2.8)
is fulfilled.

However, in practical applications, such as MRI, this is only a hypothetical
case, since the periodization $\tilde{f}$ cannot be sampled. Due to a limited
coverage of space by the acquisition, the function $f$ is typically only on a
bounded domain, w.l.o.g. for
$\boldsymbol{x}\in\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}$, cf. [1]. Thus,
we need to assume that $f$ is small outside the interval
$\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}$, such that
$\tilde{f}(\boldsymbol{x}_{j})\approx f(\boldsymbol{x}_{j})$, and we have to
deal with the approximation

$\displaystyle\hat{f}(\boldsymbol{k})\approx\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$
(3.5)

This is to say, by (3.4) the error in the approximation (3.5) solely occurs
because $f$ is not known on whole $\mathbb{R}^{d}$. Note that the
reconstruction (3.5) can also be denoted as
$\boldsymbol{\hat{f}}\approx\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$
with $\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}$ from (2.5).

######  Remark 3.1.

As known by [6, Ex. 1.22], the periodization of the function
$g(\boldsymbol{t})=\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{x}}$,
$\boldsymbol{t}\in\left[-\frac{M}{2},\frac{M}{2}\right)^{d}$, with
$\boldsymbol{x}\in\mathbb{C}^{d}$ fixed, possesses the absolutely and
uniformly convergent Fourier series

$\displaystyle\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{x}}=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{y}_{\boldsymbol{\ell}}}\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)},\quad\boldsymbol{x}\in\mathbb{C}^{d},$
(3.6)

with the $d$-variate $\mathrm{sinc}$ function
$\mathrm{sinc}(\boldsymbol{x})\coloneqq\prod_{t=1}^{d}\mathrm{sinc}(x_{t})$
with

$\displaystyle\mathrm{sinc}(x)\coloneqq\left\\{\begin{array}[]{ll}\frac{\sin
x}{x}&\quad x\in\mathbb{R}\setminus\\{0\\}\,,\\\\[4.30554pt] 1&\quad
x=0\,,\end{array}\right.$

and equispaced points
$\boldsymbol{y}_{\boldsymbol{\ell}}\coloneqq\boldsymbol{M}^{-1}\odot{\boldsymbol{\ell}}=\left(M^{-1}\ell_{1},\dots,M^{-1}\ell_{d}\right)^{T}$,
$\boldsymbol{\ell}\in\mathbb{Z}^{d}$. Since the Fourier coefficients
$\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)}$
in (3.6) are $\ell_{2}$-summable, we may introduce the $\mathrm{sinc}$
operator

$\displaystyle\mathcal{C}\coloneqq\bigg{(}\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)}\bigg{)}_{j=1,\,\boldsymbol{\ell}\in\mathbb{Z}^{d}}^{N}.$
(3.7)

By additionally defining the one-sided infinite Fourier matrix

$\displaystyle\mathcal{F}\coloneqq\left(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{y}_{\boldsymbol{\ell}}}\right)_{\boldsymbol{\ell}\in\mathbb{Z}^{d},\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}},$

the matrix product $\mathcal{C}\mathcal{F}$ can be written as

$\displaystyle\mathcal{C}\mathcal{F}=\Bigg{(}\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{y}_{\boldsymbol{\ell}}}\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)}\Bigg{)}_{j=1,\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}^{N}=\boldsymbol{A}$
(3.8)

with the nonequispaced Fourier matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in
(2.4), since the components of (3.8) coincide with point evaluations of (3.6)
at $\boldsymbol{x}=\boldsymbol{x}_{j}$, $j=1,\dots,N$, and
$\boldsymbol{t}=\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$.

###  3.2 Connection to previous work

Analogous to [4, Thm. 3.8] we aim to extend the approximation (3.4) onto the
whole interval, i. e., we consider

$\displaystyle\hat{f}(\boldsymbol{v})\approx\tilde{h}(\boldsymbol{v})\coloneqq\sum_{j=1}^{N}w_{j}\,\tilde{f}(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}},\quad\boldsymbol{v}\in\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}.$

Inserting this into the inverse Fourier transform (1.2) yields

$\displaystyle f(\boldsymbol{x})$
$\displaystyle=\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}\approx\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\tilde{h}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\sum_{j=1}^{N}w_{j}\,\tilde{f}(\boldsymbol{x}_{j})\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}(\boldsymbol{x}_{j}-\boldsymbol{x})}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\sum_{j=1}^{N}w_{j}\,\tilde{f}(\boldsymbol{x}_{j})\cdot|\mathcal{I}_{\boldsymbol{M}}|\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{x})\big{)},\quad\boldsymbol{x}\in\mathbb{R}^{d}.$

Especially, by evaluation at $\boldsymbol{x}=\boldsymbol{x}_{s}$,
$s=1,\dots,N$, we obtain

$\displaystyle f(\boldsymbol{x}_{s})$
$\displaystyle\approx\sum_{j=1}^{N}w_{j}\,\tilde{f}(\boldsymbol{x}_{j})\cdot|\mathcal{I}_{\boldsymbol{M}}|\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{s})\big{)}.$

Therefore, one could also aim to choose suitable weights $w_{j}$ based on this
equation. More precisely, by defining the nonequispaced $\mathrm{sinc}$ matrix

$\displaystyle\boldsymbol{C}_{\mathrm{n}}\coloneqq\bigg{(}\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{s})\big{)}\bigg{)}_{j,s=1}^{N}\
\in\mathbb{R}^{N\times N},$ (3.9)

as well as $\boldsymbol{f}=\left(f(\boldsymbol{x}_{j})\right)_{j=1}^{N}$ and
$\boldsymbol{\tilde{f}}=(\tilde{f}(\boldsymbol{x}_{j}))_{j=1}^{N}$, one may
try to find a weight matrix $\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}$
such that the approximation
$\boldsymbol{f}\approx|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}\boldsymbol{\tilde{f}}$
is best as possible. Note that still $\boldsymbol{\tilde{f}}$ is unknown, such
that we have to deal with the overall approximation
$\boldsymbol{f}\approx|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}\boldsymbol{f}$,
i. e., one would ideally aim for
$\boldsymbol{I}_{N}=|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}$.
We remark that on the main diagonal
$|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}=\boldsymbol{I}_{N}$
reads as

$\displaystyle\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}}|}=w_{j}\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{j})\big{)}=w_{j}\,\mathrm{sinc}(0)=w_{j},\quad
j=1,\dots,N.$

However, for all other entries with $j\neq s$ we would need
$\mathrm{sinc}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{s}))=0$, which is only
true for $\boldsymbol{x}_{j}$ equispaced. In other words, for arbitrary points
$\boldsymbol{x}_{j}$ equality in
$|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}=\boldsymbol{I}_{N}$
cannot be fulfilled for any weights.

Hence, one can only look for an approximate solution, e. g. by considering the
least squares problem

$\displaystyle\underset{\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}}{\text{Minimize
}}\
\||\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}-\boldsymbol{I}_{N}\|_{\mathrm{F}}^{2}.$

When denoting the $j$-th column of $\boldsymbol{W}$, $\boldsymbol{I}_{N}$ and
$\boldsymbol{C}_{\mathrm{n}}$ as $\boldsymbol{w}_{j}$, $\boldsymbol{e}_{j}$
and $[\boldsymbol{C}_{\mathrm{n}}]_{j}$, respectively, we may rewrite the
Frobenius norm by only considering the nonzero entries via

$\displaystyle\||\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{W}-\boldsymbol{I}_{N}\|_{\mathrm{F}}^{2}$
$\displaystyle=\sum_{j=1}^{N}\||\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{C}_{\mathrm{n}}\boldsymbol{w}_{j}-\boldsymbol{e}_{j}\|_{2}^{2}=\sum_{j=1}^{N}\left\||\mathcal{I}_{\boldsymbol{M}}|\cdot[\boldsymbol{C}_{\mathrm{n}}]_{j}\,{w_{j}}-\boldsymbol{e}_{j}\right\|_{2}^{2}.$

Thus, as stated in [2] (without proof) the least squares solution to the
minimization problem is given by

$\displaystyle{w_{j}}$
$\displaystyle=\frac{1}{|\mathcal{I}_{\boldsymbol{M}}|}\,[\boldsymbol{C}_{\mathrm{n}}]_{j}^{\dagger}\,\boldsymbol{e}_{j}=\frac{1}{|\mathcal{I}_{\boldsymbol{M}}|}\left([\boldsymbol{C}_{\mathrm{n}}]_{j}^{*}\,[\boldsymbol{C}_{\mathrm{n}}]_{j}\right)^{-1}[\boldsymbol{C}_{\mathrm{n}}]_{j}^{*}\,\boldsymbol{e}_{j}$
$\displaystyle=\frac{\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{j})\big{)}}{|\mathcal{I}_{\boldsymbol{M}}|\,\sum_{s=1}^{N}\,\mathrm{sinc}^{2}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{s}))}=\frac{1}{|\mathcal{I}_{\boldsymbol{M}}|}\left(\sum_{s=1}^{N}\mathrm{sinc}^{2}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{s}))\right)^{-1}.$
(3.10)

######  Remark 3.2.

Note that in [2] it was claimed that this approach coincides with the one in
[7] only considering finite sections of (3.7) with
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$. However, we remark that
this claim only holds asymptotically for
$|\mathcal{I}_{\boldsymbol{M}}|\to\infty$.

By the classical sampling theorem of Shannon-Whittaker-Kotelnikov, see [10, 9,
5], any bandlimited function $f\in L_{2}(\mathbb{R}^{d})$ with maximum
bandwidth $\boldsymbol{M}$ can be recovered from its uniform samples
$f(\boldsymbol{y}_{\boldsymbol{\ell}})$, $\boldsymbol{\ell}\in\mathbb{Z}^{d}$,
and we have

$\displaystyle
f(\boldsymbol{x})=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}f(\boldsymbol{y}_{\boldsymbol{\ell}})\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)},\quad\boldsymbol{x}\in\mathbb{R}^{d}.$
(3.11)

Now we apply the sampling theorem of Shannon-Whittaker-Kotelnikov to the
shifted $\mathrm{sinc}$ function
$f(\boldsymbol{x})=\mathrm{sinc}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}))\in
L_{2}(\mathbb{R}^{d})\cap C(\mathbb{R}^{d})$ with $j$ fixed. By evaluation at
$\boldsymbol{x}=\boldsymbol{x}_{s}$, $s=1,\dots,N$, we obtain

$\displaystyle\mathrm{sinc}\big{(}$ $\displaystyle
M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}_{s})\big{)}=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)}\,\mathrm{sinc}\big{(}M\pi(\boldsymbol{x}_{s}-\boldsymbol{y}_{\boldsymbol{\ell}})\big{)}.$
(3.12)

Using the $\mathrm{sinc}$ operator $\mathcal{C}$ in (3.7) and the
nonequispaced $\mathrm{sinc}$ matrix in (3.9), this can be written as

$\displaystyle\boldsymbol{C}_{\mathrm{n}}=\mathcal{C}\mathcal{C}^{*}\in\mathbb{R}^{N\times
N}.$ (3.13)

Hence, a restriction to finitely many $\boldsymbol{\ell}$ in (3.12)
corresponds to uniform truncation of a Shannon series (3.11), which is known
as a poor approximation due to the slow convergence of the $\mathrm{sinc}$
function. Thus, equality in (3.13) is only satisfied for the operator
$\mathcal{C}$ in (3.7), while considering finite sections of (3.7) with
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$ implies a poor
approximation of $\boldsymbol{C}_{\mathrm{n}}$.

##  4 Numerics

(a) Spiral grid

(b) Polar (blue) and modified polar (red) grid

(c) Logarithmic modified polar grid

Figure 4.1:  Exemplary grids of size $R=12$ and $T=2R$.

######  Example 4.1.

Analogous to [4, Ex. 5.4] we firstly consider the reconstruction of
trigonometric polynomials using the example of the Shepp-Logan phantom, see
Fig. 2(a), i. e., the phantom data shall be treated as Fourier coefficients
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
of a trigonometric polynomial (1.1). For given
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{2}$,
$j=1,\dots,N$, we then compute the evaluations $f(\boldsymbol{x}_{j})$ of
(1.1) by means of an NFFT and use the resulting vector as input for the
reconstruction.

Since for $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$ the optimality of weights
computed by (2.6) was shown in [4], we now consider
$|\mathcal{I}_{\boldsymbol{2M}}|>N$ and compare the density compensation
factors mentioned in Section 2. More precisely, for the spiral grid, cf. Fig.
1(a), of size $R=M$, $T=2R$, we study (2.6), (2.7) and (2.9). The resulting
reconstruction of the phantom of size $M=64$ is presented in Fig. 4.2 (top)
including a detailed view of the 52nd row (bottom). It can easily be seen that
since the exactness condition $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$ (see
Theorem 2.1) is violated, the weights computed by (2.6) do not yield an exact
reconstruction, cf. Fig. 2(b). Note that the results using the least squares
approximation by (2.7) are just as bad, cf. Fig. 2(c). Merely, for the weights
computed by (2.9) we see a slight improvement in Fig. 2(d), since there are no
artifacts on the inside of the phantom and the differences to the original
phantom are not visible anymore. We remark that the comparatively small choice
of $M=64$ is necessary for (2.9) being computationally affordable.

(a) Phantom

(b) Use of (2.6)

(c) Use of (2.7)

(d) Use of (2.9)

Figure 4.2:  Reconstruction of the Shepp-Logan phantom of size $M=64$ (top)
via density compensation factors computed by (2.6), (2.7) and (2.9) for the
spiral grid, cf. Fig. 1(a), of size $R=M$, $T=2R$; as well as a detailed view
of the 52nd row (bottom).

######  Example 4.2.

Summing up, we examine the reconstruction properties for bandlimited functions
$f\in L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ with maximum bandwidth
$\boldsymbol{M}$. To determine the errors properly, we firstly specify a
compactly supported function $\hat{f}$ and consequently compute its inverse
Fourier transform (1.2), such that its samples $f(\boldsymbol{x}_{j})$ for
given $\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{2}$,
$j=1,\dots,N$, can be used for the reconstruction of the samples
$\hat{f}(\boldsymbol{k})$, $\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$. As
in [4, Ex. 5.5] we consider the tensorized function
$\hat{f}(\boldsymbol{v})=g(v_{1})\cdot g(v_{2})$, where $g(v)$ is the one-
dimensional triangular pulse
$g(v)\coloneqq(1-\left|\tfrac{v}{b}\right|)\cdot\chi_{[-b,b]}(v)$. Then for
all $b\in\mathbb{N}$ with $b\leq\frac{M}{2}$ the associated inverse Fourier
transform

$\displaystyle
f(\boldsymbol{x})=\int_{\mathbb{R}^{2}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}=b^{2}\,\mathrm{sinc}^{2}(b\pi\boldsymbol{x}),\
\;\boldsymbol{x}\in\mathbb{R}^{2},$ (4.1)

is bandlimited with bandwidth $\boldsymbol{M}$. For this test function, we
compare the density compensation methods from Section 3, i. e., we study the
computaion schemes (2.6), (2.9) and (3.2).

As a first experiment we fix $M=32$ and $b=12$ and consider the case
$|\mathcal{I}_{\boldsymbol{2M}}|\leq N$, which yields optimality for
trigonometric polynomials. In addition to the real-world sampling of (4.1)
also examine the artificial sampling data

$\displaystyle\tilde{f}(\boldsymbol{x}_{j})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}({\boldsymbol{k}})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}$
(4.2)

of the periodization (3.2). A visualization of the chosen test function $f$
and the difference $f-\tilde{f}$ can be found in Fig. 4.3.

(a) Test function $f$

(b) $f-\tilde{f}$

Figure 4.3:  The test function (4.1) and its periodization (4.2).

For the modified polar grid, cf. Fig. 1(b), of size $R=2M$, $T=2R$, we use
these two kinds of sampling data to compute the reconstructions
$\boldsymbol{\tilde{h}}\coloneqq({\tilde{h}}^{\mathrm{w}}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
cf. (2.2), and the pointwise errors
$\big{|}\boldsymbol{\tilde{h}}-\boldsymbol{\hat{f}}\big{|}$. The corresponding
results are displayed in Fig. 4.4. It can easily be seen that for the
artificial sampling data (4.2) the weights computed by (2.6) indeed yield an
exact reconstruction, see Fig. 4(a) (bottom), and thus are optimal. However,
in the more realistic setting the results are not as good, but the weights by
(2.6) and (2.9), see Fig. 4(a) and 4(b) (top), produce nearly the same error
as a reconstruction on an equispaced grid. Hence, for bandlimited functions
the truncation $\tilde{f}(\boldsymbol{x}_{j})\approx f(\boldsymbol{x}_{j})$ is
the dominating error term and therefore reconstruction errors smaller than the
ones shown in Fig. 4.4 cannot be expected. In addition, we remark that (3.2)
is not as good, for none of the sampling data, see Fig. 4(c). Note that we
chose comparatively small $M=32$ in order that the computation schemes (2.9)
and (3.2) are affordable.

(a) Use of (2.6)

(b) Use of (2.9)

(c) Use of (3.2)

Figure 4.4:  Pointwise error of the reconstruction of the tensorized
triangular pulse with $M=32$ and $b=12$, via density compensation factors
computed by (2.6), (2.9) and (3.2) for the modified polar grid, cf. Fig. 1(b),
of size $R=2M$, $T=2R$, using samples $f(\boldsymbol{x}_{j})$ (top) and
artificial samples $\tilde{f}(\boldsymbol{x}_{j})$ (bottom).

In a second experiment we sample the function (4.1) with $M=64$ and $b=24$ at
logarithmic modified polar grids, cf. Fig. 1(c), of different sizes $R$,
$T=2R$. The corresponding relative errors

$\displaystyle\frac{\|\boldsymbol{\tilde{h}}-\boldsymbol{\hat{f}}\|_{2}}{\|\boldsymbol{\hat{f}}\|_{2}}$
(4.3)

can be found in Table 4.1. Note that for $M=64$ we have
$|\mathcal{I}_{\boldsymbol{M}}|=4096$ and
$|\mathcal{I}_{\boldsymbol{2M}}|=16384$. We observe that for
$|\mathcal{I}_{\boldsymbol{2M}}|\leq N$ the errors from (2.6) and (2.9) are
comparably good, while except for $|\mathcal{I}_{\boldsymbol{M}}|>N$ the
scheme (2.9) produces optimal results.

R | N |  Use of (2.6) |  Use of (2.9) |  Use of (3.2)
---|---|---|---|---
40 | 3565 | 4.4908e-01 | 1.7608e-01 | 2.0475e-01
48 | 5145 | 1.0886e-01 | 2.0690e-02 | 1.5829e-01
56 | 7149 | 3.6632e-02 | 8.0215e-03 | 1.5401e-01
64 | 9429 | 2.5109e-02 | 4.7988e-03 | 1.8337e-01
72 | 11965 | 7.6871e-03 | 4.1096e-03 | 2.0633e-01
80 | 14909 | 5.5991e-03 | 3.8507e-03 | 2.1932e-01
88 | 18153 | 3.8889e-03 | 3.9853e-03 | 2.2665e-01
96 | 21589 | 4.2240e-03 | 3.7917e-03 | 2.3092e-01
Table 4.1:  Relative errors (4.3) of the reconstruction of the tensorized
triangular pulse with $M=64$ and $b=24$, via density compensation factors
computed by (2.6), (2.9) and (3.2) for log. modified polar grids, cf. Fig.
1(c), of different sizes.

## Acknowledgments

Melanie Kircheis acknowledges the support from the BMBF grant 01$\mid$S20053A
(project SA$\ell$E). Daniel Potts acknowledges the funding by Deutsche
Forschungsgemeinschaft (German Research Foundation) – Project–ID 416228727 –
SFB 1410.

## References

  * [1] H. Eggers, M. Kircheis, and D. Potts.  Non-Cartesian MRI reconstruction.  In M. Doneva, M. Akcakaya, and C. Prieto, editors, Magnetic Resonance Image Reconstruction: Theory, Methods and Applications, volume 6. Academic Press, 2022.
  * [2] L. Greengard, J.-Y. Lee, and S. Inati.  The fast sinc transform and image reconstruction from nonuniform samples in $k$-space.  Commun. Appl. Math. Comput. Sci., 1:121–131, 2006.
  * [3] D. Huybrechs.  Stable high-order quadrature rules with equidistant points.  J. Comput. Appl. Math., 231(2):933–947, 2009.
  * [4] M. Kircheis and D. Potts.  Fast and direct inversion methods for the multivariate nonequispaced fast Fourier transform.  arXiv: 2302.00414, 2023.
  * [5] V. A. Kotelnikov.  On the transmission capacity of the “ether” and wire in electrocommunications.  In Modern Sampling Theory: Mathematics and Application, pages 27–45. Birkhäuser, Boston, 2001.  Translated from Russian.
  * [6] J. Lund and K. L. Bowers.  Sinc Methods for Quadrature and Differential Equations.  Society for Industrial and Applied Mathematics, 1992.
  * [7] J. G. Pipe and P. Menon.  Sampling density compensation in MRI: rationale and an iterative numerical solution.  Magn. Reson. Med., 41:179–186, 1999.
  * [8] G. Plonka, D. Potts, G. Steidl, and M. Tasche.  Numerical Fourier Analysis.  Applied and Numerical Harmonic Analysis. Birkhäuser, 2018.
  * [9] C. E. Shannon.  Communication in the presence of noise.  Proc. I.R.E., 37:10–21, 1949.
  * [10] E. T. Whittaker.  On the functions which are represented by the expansions of the interpolation theory.  Proc. R. Soc. Edinb., 35:181–194, 1915.


# Network Multibaker Maps, Information Theory and Stochastic Thermodynamics

Bernhard Altaner  Jürgen Vollmer  Max-Planck-Institut für Dynamik und
Selbstorganisation (MPIDS), 37077 Göttingen, Germany  Fakultät für Physik,
Universität Göttingen, 37077 Göttingen, Germany

###### Abstract

Stochastic thermodynamics (ST) summarizes recent advances in nonequilibrium
statistical mechanics of systems featuring local thermodynamic equilibrium.
However, the mathematical description of ST is largely independent of this
thermodynamic context. Here, we introduce novel network multibaker maps as a
model system for reversible, chaotic dynamics. As a result, the mathematical
formalism of ST emerges from information theory, and is shown to be consistent
with previous work.

###### pacs:

05.45.-a, 05.70.Ln,05.40.-a, 89.70.Cf

## Introduction

_Stochastic thermodynamics_ (ST) is a recent approach to the statistical
mechanics of nonequilibrium systems Seifert (2008). Using both random
variables and ensemble properties it extends the notions of entropy and
entropy production to stochastic trajectories Maes (2004); Seifert (2005).
Using the notion of time-reversal, ST provides elegant derivations Seifert
(2012) of stochastic fluctuation Kurchan (1998); Lebowitz and Spohn (1999) and
work relations Jarzynski (1997); Crooks (1999). It also yields universal
results on thermodynamic efficiency of small machines Esposito  _et al._
(2009) and the conversion of information to work Horowitz and Parrondo (2011);
Sagawa and Ueda (2012). Moreover, advanced experimental and simulation
techniques allow direct testing of these theoretical results Toyabe  _et al._
(2010); Tietz  _et al._ (2006); Ciliberto  _et al._ (2010). ST relies on a
thermodynamically consistent interpretation of the stochastic process within a
framework of an assumed thermodynamic _local equilibrium_. Surprisingly, the
mathematical formulation of ST seems not to depend on the specific physical
context Seifert (2011); Esposito (2012). What is needed is the _Markov
property_ of the dynamics, though for certain relations even this assumptions
can be relaxed Speck and Seifert (2007).

In this work, we introduce the notion of _network multibaker maps_ (NMBM) and
take a dynamical system’s perspective on ST. Multibaker maps serve as a
paradigm for studies on the relation of microscopic dynamics, statistical
physics and transport phenomena Vollmer (2002); Colangeli  _et al._ (2011).
They fulfill many of the generic features of chaotic dynamics: Reversibility,
hyberbolicity, homo- and heteroclinic orbits, a natural Markov partition and
ergodicity Hopf (1948); Gaspard (2005). Yet, their simple structure admits an
explicit calculation of the evolution of ensembles which are represented by
non-stationary phase-space densities Vollmer  _et al._ (1997); Vollmer (2002).
Here, we explicitly calculate the full evolution of information-theoretical
quantities for NMBM. We distinguish three fundamental contributions, which are
in one-to-one correspondance to the entropy expressions used in ST. Further,
we show that our approach is consistent with previous work based on phase-
space contraction Rondoni  _et al._ (2000); Evans and Searles (2002).

## Network multibaker maps

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a connected graph with $N$
vertices $i\in\mathcal{V}:=[1,2,\ldots,N]$ and edges
$e\in\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$. Vertices represent the
fundamental states of some effective (for instance, thermodynamic)
description. Edges indicate possible transitions between these states. We
assume bi-directional edges, i.e.
$e:=(i,j)\in\mathcal{E}\Leftrightarrow\overline{e}:=(j,i)\in\mathcal{E}$, and
that each vertex is connected to itself, i.e.  $(i,i)\in\mathcal{E},\,\forall
i\in\mathcal{V}$. Let
$\mathcal{V}_{i}:={\\{j\in\mathcal{V}:(i,j)\in\mathcal{E}\\}}$ be the set of
vertices that are connected to a state $i$. Denote by
$k_{i}:=\left|\mathcal{V}_{i}\right|$ the degree of vertex $i$. Associate with
each vertex $i$ a rectangular _cell_
$\mathcal{C}_{i}:=[0,1]\times\Pi_{i}\cdot[0,1]$ with area (Lebesgue measure)
$\Pi_{i}$ consisting of the set of phase-space points $x\in\mathcal{C}_{i}$
which belong to state $i$. The overall _phase space_ $\Gamma$ of the system is
the disjoint union $\Gamma:=\cup_{i=1}^{N}\mathcal{C}_{i}$.

Figure 1: Dynamics of the network multibaker map. For each of the $k_{i}$
neighbours vertices $j\in\mathcal{V}_{j}$ of vertex $i$, a horizontal strip
$C^{i}_{j}\subset\mathcal{C}_{i}$ of relative height $s^{i}_{j}$ is affine-
linearly mapped to a vertical strip
$\tilde{}\mathcal{C}^{i}_{j}\subset\mathcal{C}_{j}$ of width
$\tilde{s}^{i}_{j}$.

A NMBM, $\Phi:\Gamma\to\Gamma$, deterministically maps phase-space points to
adjacent cells. It is specified geometrically (cf. Fig. 1) by dividing each
cell $C_{i}$ into $k_{i}$ horizontal strips of finite relative height
$s^{i}_{j}>0$. More precisely, with the offset
$b^{i}_{j}=\sum_{k<j}s^{i}_{k}$:

$\mathcal{C}^{i}_{j}:=\begin{cases}[0,1]\times\Pi_{i}\cdot[b^{i}_{j},b^{i}_{j}+s^{i}_{j})&i\in\mathcal{V},j\in\mathcal{V}_{i},\\\
\emptyset&\text{else}.\end{cases}$ (1)

The dynamics maps each horizontal strip
$\mathcal{C}^{i}_{j}\subset\mathcal{C}_{i}$ to a vertical strip
$\tilde{}\mathcal{C}^{i}_{j}:=\Phi\mathcal{C}^{i}_{j}\subset\mathcal{C}_{j}$
in an affine-linear way:

$\tilde{}\mathcal{C}^{i}_{j}:=[\tilde{b}^{i}_{j},\tilde{b}^{i}_{j}+\tilde{s}^{i}_{j})\times\Pi_{j}\cdot[0,1]$
(2)

Analogously, $\tilde{s}^{i}_{j}$ denote the relative width of the vertical
strips and the offsets read $\tilde{b}^{i}_{j}=\sum_{k<i}\tilde{s}^{i}_{j}$.
Points in any strip are mapped affine-linearly, such that the horizontal
direction is contracted by a factor $\tilde{s}^{i}_{j}<1$ whereas the vertical
direction is expanded by a factor $(s^{i}_{j})^{-1}>1$. The multibaker map is
fully defined by specifying $\mathcal{G}$ and the numbers $\Pi_{i}$,
$s^{i}_{j}$ and $\tilde{s}^{i}_{j}$.

In deriving fluctuation relations in both stochastic and deterministic
dynamics, one applies a natural notion of time-reversal Lebowitz and Spohn
(1999); Rondoni  _et al._ (2000); Maes and Netočný (2003); Seifert (2005)
which has its origins in the time-reversibility of fundamental physical laws.
Following Ref. Tel and Vollmer (2000) we are interested in dynamics with a
transformation $\mathcal{I}$ on $\Gamma$ which a) is an involution, b)
preserves phase-space volume $\left|\cdot\right|$ and c) facilitates time-
reversal:

$\displaystyle\mathcal{I}\circ\mathcal{I}$ $\displaystyle=\text{id},$ (3a)
$\displaystyle\left|\mathcal{I}^{-1}A\right|$
$\displaystyle=\left|A\right|,~{}\forall A\subset\Gamma,$ (3b)
$\displaystyle\Phi^{-1}x$
$\displaystyle=\left(\mathcal{I}\circ\Phi\circ\mathcal{I}\right)x,~{}\forall
x\in\Gamma.$ (3c)

_Properly thermostatted multibaker maps_ Tel and Vollmer (2000) obey

$\displaystyle\tilde{s}^{i}_{j}=s^{j}_{i},\,\forall i,j,$ (4)

such that Eq. (3) is fulfilled by the involution Vollmer (2002)

$\displaystyle\mathcal{I}:(x,y)$
$\displaystyle\mapsto(1-\Pi_{i}^{-1}y,\Pi_{i}(1-x)),\,\forall(x,y)\in\mathcal{C}_{i},$
(5)

which additionally is local each cell (cf.  Fig. 2), i.e.

$\mathcal{I}\mathcal{C}_{i}=\mathcal{C}_{i},\,\forall i$ (6)

With the general definition
$\tilde{}\mathcal{C}^{i}_{j}:=\Phi\mathcal{C}_{i}\cap\mathcal{C}_{j}$, Eq. (4)
follows solely from Eqs. (3) and (6), which ensure that
$\tilde{}\mathcal{C}^{i}_{j}$ and $\mathcal{C}^{j}_{i}$ are interchanged by
the involution (cf.  Fig. 2).

Figure 2:  The multibaker dynamics feature a volume preserving involution
$\mathcal{I}$, cf.  Eq. (5), which interchanges horizontal and vertical
stripes in one cell, i.e.
$\mathcal{I}\tilde{}\mathcal{C}^{i}_{j}=\mathcal{C}^{j}_{i}$ The dynamics
$\Phi$ is reversible, but generally not volume preserving.

## Densities and trajectories

In statistical physics, the probability of finding a system at some point
$x\in\Gamma$ is given by a time-dependent ensemble measure. We assume that
this ensemble is absolutely continuous with respect to the Lebesgue measure,
i.e. it is characterized by a phase-space density
$\varrho^{\left(\nu\right)}(x)$, where the superscript (ν) indicates discrete
time. This phase-space density evolves according to $\Phi$ and generically
develops an intricate structure as the number of iterations increases.
Practically, preparation and resolution of phase-space densities is restricted
to a coarse-grained level, which is represented by the vertices of
$\mathcal{G}$. Hence, we define a coarse-grained density
$\rho^{\left(\nu\right)}(x)$ to take the average value
$\rho^{\left(\nu\right)}_{i}$ of $\varrho^{\left(\nu\right)}$ on cell
$\mathcal{C}_{i}$:

$\displaystyle\rho_{i}^{\left(\nu\right)}$
$\displaystyle:=\Pi_{i}^{-1}\int_{\mathcal{C}_{i}}\varrho^{\left(\nu\right)}(x)\,\mathrm{d}x,$
(7a) $\displaystyle\rho^{\left(\nu\right)}(x)$
$\displaystyle:=\sum_{i}\chi_{i}(x)\rho^{\left(\nu\right)}_{i},$ (7b)

where $\chi_{i}$ is the characteristic or indicator function of
$\mathcal{C}_{i}$.

Due to the preparation procedure, fine- and coarse-grained ensembles agree
initially, i.e. $\varrho^{(0)}(x)\equiv\rho^{(0)}(x)$. To follow the evolution
of $\varrho^{(0)}$ to $\varrho^{(\nu)}$ we introduce trajectories
$\underline{\omega}$ as $(\nu+1)$-tuples with entries
$\omega_{k}\in\mathcal{V}$:

$\underline{\omega}:=\left(\omega_{0},\omega_{1},\ldots,\omega_{\nu}\right).$
(8)

The entry $\omega_{k}$ denotes the vertex visited at time $k$. Hence, a
trajectory is a finite string of a the symbolic sequence associated with the
(natural) Markov partition Beck and Schögl (1995).

Let the set $\tilde{}\mathcal{C}\left[\underline{\omega}\right]$ denote all
phase-space points that move together along a fixed trajectory
$\underline{\omega}$ and hence reside in cell $\mathcal{C}_{\omega_{\nu}}$ at
time $\nu$. Denoting by $\underline{\omega}^{(k)}$ the first $k\leq\nu$ steps
of $\underline{\omega}$, we define recursively:

$\displaystyle\tilde{}\mathcal{C}[\underline{\omega}^{(k)}]:=\Phi(\tilde{}\mathcal{C}[\underline{\omega}^{(k-1)}]\cap\mathcal{C}^{\omega_{k-1}}_{\omega_{k}}),\quad\forall\,0<k\leq\nu,$
(9)

where
$\tilde{}\mathcal{C}\left[\underline{\omega}^{(0)}\right]:=\mathcal{C}_{\omega_{0}}$.
One can easily verify that $\tilde{}\mathcal{C}[\underline{\omega}^{(k-1)}]$
is a vertical strip. Intersecting with the horizontal strip
$\mathcal{C}^{\omega_{k-1}}_{\omega_{k}}$ reduces its volume by a factor
$s^{\omega_{k-1}}_{\omega_{k}}<1$. Through contraction and expansion in
horizontal and vertical directions, respectively, $\Phi$ changes the volume by
another factor
$\tilde{s}^{\omega_{k-1}}_{\omega_{k}}/s^{\omega_{k-1}}_{\omega_{k}}$. Hence,
$\left|\mathcal{C}\left[\underline{\omega}^{(k)}\right]\right|=\tilde{s}^{\omega_{k-1}}_{\omega_{k}}\left|\mathcal{C}\left[\underline{\omega}^{(k-1)}\right]\right|$
and after iteration to $k=\nu$,

$\displaystyle\left|\mathcal{C}[\underline{\omega}^{(\nu)}]\right|=\Pi_{\omega_{\nu}}\prod_{k=1}^{\nu}\tilde{s}^{\omega_{k-1}}_{\omega_{k}}=\Pi_{\omega_{\nu}}\prod_{k=1}^{\nu}s^{\omega_{k}}_{\omega_{k-1}},$
(10)

where for the last equality we used Eq. (4).

Because the initial microscopic density $\varrho^{(0)}$ is uniform on each
cell, so is the density on $\tilde{}\mathcal{C}[\underline{\omega}]$ for any
$\underline{\omega}$. Denote by $\varrho[\underline{\omega};k]$ the density on
the $(\nu-k)$th pre-image
$\Phi^{\nu-k}\tilde{}\mathcal{C}[\underline{\omega}]$. Probability
conservation and uniform contraction imply

$\displaystyle\varrho[\underline{\omega};k]=(\eta^{\omega_{k-1}}_{\omega_{k}})^{-1}\varrho[\underline{\omega};k-1]$
(11)

where we introduced the contraction factor

$\eta^{\omega_{k-1}}_{\omega_{k}}:=\frac{\left|\tilde{}\mathcal{C}^{\omega_{k-1}}_{\omega_{k}}\right|}{\left|\mathcal{C}^{\omega_{k-1}}_{\omega_{k}}\right|}\equiv\frac{\Pi_{\omega_{k}}\tilde{s}^{\omega_{k-1}}_{\omega_{k}}}{\Pi_{\omega_{k-1}}s^{\omega_{k-1}}_{\omega_{k}}}.$
(12)

Iterating this condition, we find with Eq. (4) that the density
$\varrho^{(\nu)}(x)$ at a point $x$ with history $\underline{\omega}(x)$
obeys:

$\varrho^{(\nu)}(x)=\varrho\left[\underline{\omega}(x);\nu\right]=\rho_{\omega_{0}}\prod_{k=1}^{\nu}\left[\frac{s^{\omega_{k-1}}_{\omega_{k}}}{s^{\omega_{k}}_{\omega_{k-1}}}\right]\frac{\Pi_{\omega_{0}}}{\Pi_{\omega_{\nu}}}.$
(13)

Integrating $\varrho(x)$ over a cell $\Pi_{i}$ yields the probability

$\displaystyle
p_{i}^{(k)}:=\int_{\mathcal{C}_{i}}\varrho^{\left(k\right)}(x)\,\,\mathrm{d}x\,=\Pi_{i}\rho_{i}^{(k)},$
(14)

which evolves according to the _Master equation_ Beck and Schögl (1995)

$\displaystyle p_{i}^{(k)}=\sum_{j}\left[p_{j}^{(k-1)}s^{j}_{i}\right].$ (15)

This establishes the connection to _Markov chains_ and hence the (discrete)
formulation of ST Seifert (2012).

## Information and entropy

The Shannon entropy Shannon (1948) of a probability density $\varrho(x)$ on
$\Gamma$ is defined as

$\displaystyle\mathcal{S}\left[\varrho\right]:=-\int_{\Gamma}\varrho\,\log\varrho\,\,\mathrm{d}x.$
(16)

With the results (10) and (13) we can calculate the fine-grained entropy of
$\varrho^{\left(\nu\right)}(x)$ as a sum over all trajectories
$\underline{\omega}^{\left(\nu\right)}$ of length $\nu$ weighted with
$\left|\mathcal{C}[\underline{\omega}^{\left(\nu\right)}]\right|$

$\displaystyle
S_{\mathrm{fg}}^{(\nu)}:=\mathcal{S}\left[\varrho^{(\nu)}\right]=-\sum_{\underline{\omega}^{\left(\nu\right)}}\left[\left|\mathcal{C}[\underline{\omega}^{\left(\nu\right)}]\right|\times\left(\varrho[\underline{\omega};\nu]\log\varrho[\underline{\omega};\nu]\right)\right]$
$\displaystyle=\sum_{\underline{\omega}^{\left(\nu\right)}}\left[\mathbb{P}[\underline{\omega}]\left(\log\Pi_{\omega_{\nu}}-\log
p^{(0)}_{\omega_{0}}-\sum_{k=1}^{\nu}B^{\omega_{k-1}}_{\omega_{k}}\right)\right]$
(17)

where

$\displaystyle\mathbb{P}[\underline{\omega}^{\left(\nu\right)}]$
$\displaystyle\equiv\left|\mathcal{C}[\underline{\omega}^{\left(\nu\right)}]\right|\varrho[\underline{\omega};\nu]=p^{\left(0\right)}_{\omega_{0}}\prod_{k=1}^{\nu}s^{\omega_{k-1}}_{\omega_{k}},$
(18) $\displaystyle B^{i}_{j}$
$\displaystyle:=\log\,\left({s^{i}_{j}}/{s^{j}_{i}}\right).$ (19)

Similarly, we find

$\displaystyle S_{\mathrm{cg}}^{(\nu)}$
$\displaystyle:=\mathcal{S}\left[\varrho^{(\nu)}\right]=\sum_{\underline{\omega}^{\left(\nu\right)}}\left[\mathbb{P}[\underline{\omega}^{\left(\nu\right)}]\left(\log\Pi_{\omega_{\nu}}-\log
p^{(\nu)}_{\omega_{\nu}}\right)\right].$ (20)

## Trajectory functionals

The forms of Eqs. (17) and (20) suggest to write them as weighted averages of
trajectory functionals $s[\underline{\omega};\nu]$,i.e. :

$\displaystyle\left\langle\\!\left\langle
s\right\rangle\\!\right\rangle^{(\nu)}=\sum_{\underline{\omega}^{\left(\nu\right)}}\left[\mathbb{P}[\underline{\omega}]\,s[\underline{\omega};\nu]\right]$
(21)

With that we can write $S_{\mathrm{fg}}^{(\nu)}=\left\langle\\!\left\langle
s_{\mathrm{fg}}\right\rangle\\!\right\rangle^{\left(\nu\right)}$ and
$S_{\mathrm{cg}}^{(\nu)}=\left\langle\\!\left\langle
s_{\mathrm{cg}}\right\rangle\\!\right\rangle^{\left(\nu\right)}$, where the
functionals

$\displaystyle s_{\mathrm{fg}}[\underline{\omega},\nu]$
$\displaystyle=s_{\mathrm{int}}[\underline{\omega},\nu]+s_{\mathrm{vis}}[\underline{\omega},0]-s_{\mathrm{mot}}[\underline{\omega},\nu],$
(22a) $\displaystyle s_{\mathrm{cg}}[\underline{\omega},\nu]$
$\displaystyle=s_{\mathrm{int}}[\underline{\omega},\nu]+s_{\mathrm{vis}}[\underline{\omega},\nu],$
(22b)

are linear combinations of three fundamental functionals:

$\displaystyle s_{\mathrm{int}}[\underline{\omega},\nu]$
$\displaystyle=s_{\mathrm{int}}(\omega_{\nu})=\log\Pi_{\omega_{\nu}},$ (23a)
$\displaystyle s_{\mathrm{vis}}[\underline{\omega},\nu]$
$\displaystyle=s_{\mathrm{vis}}(\omega_{\nu},\nu)=-\log
p^{\left(\nu\right)}_{\omega_{\nu}},$ (23b) $\displaystyle
s_{\mathrm{mot}}[\underline{\omega},\nu]$
$\displaystyle=s_{\mathrm{mot}}[\underline{\omega}]=\sum_{k=1}^{\nu}B^{\omega_{k-1}}_{\omega_{k}}.$
(23c)

Note that the functionals $s_{\mathrm{int}}$ and $s_{\mathrm{vis}}$ depend on
the final state $\omega_{\nu}$ of the trajectory only, i.e.
$s[\underline{\omega};\nu]=s(\omega_{\nu};\nu)$. Hence, we say they have a
_local form_. Trajectory averages of local forms reduce to ensemble averages
via Eq. (15):

$\displaystyle\left\langle\\!\left\langle
s\right\rangle\\!\right\rangle^{(\nu)}=\left\langle
s\right\rangle^{\left(\nu\right)}:=\sum_{i}\left[p_{i}^{\left(\nu\right)}s(i;\nu)\right].$
(24)

Henceforth, we regard NMBM as a model for thermostatted reversible dynamics
Jepps and Rondoni (2010). The Markovian description (15) for the coarse-
grained ensemble on the discrete state establishes the connection to ST.

The _intrinsic entropy_ $s_{\mathrm{int}}$ accounts for the entropy of the
unobservable microscopic configurations. In the present dynamical context it
is calculated à la Boltzmann as the logarithm of a phase-space volume (23a).
In the thermodynamic context, it is the entropy of the local equilibrium
ensemble which reflects the state of the environment Seifert (2011). Hence,
the intrinsic entropy can be considered a “state variable” and in this respect
is closest to standard thermodynamics.

The entropy that is operationally accesible in experiments is the _visible
entropy_ $s_{\mathrm{vis}}$. It only depends on the present ensemble, i.e.
the $p_{i}^{(\nu)}$, which can be sampled by performing many measurements on
identical systems. Though it has local form and averages can be calculated
with Eq. (24), it depends on the ensemble and hence is not a state variable.

The conceptually most difficult term is the third functional,
$s_{\mathrm{mot}}$ related to the quantity $B^{i}_{j}$. In analogy to the
electromotance for electrical circuits Altaner  _et al._ (2012), we call
$B^{i}_{j}$ the _motance_ of a transition $i\to j$. Physically, a finite
motance arises from a system-environment interaction. Unlike the affinity
$A^{i}_{j}=\log({p_{i}s^{i}_{j}}/({p_{j}s^{j}_{i}}))$, the motance $B^{i}_{j}$
only depends on the dynamics and not on the ensemble Schnakenberg (1976).
Consequently, we call $s_{\mathrm{mot}}$ the _accumulated motance_.

## Entropy variation and production

For each fundamental functional $s[\underline{\omega},\nu]$ there is an
associated fundamental variation
$\sigma[\underline{\omega},\nu]:=s[\underline{\omega},\nu]-s[\underline{\omega}^{(\nu-1)},\nu-1]$.
The fundamental variations

$\displaystyle\sigma_{\mathrm{int}}[\underline{\omega},\nu]$
$\displaystyle=\log{\Pi_{\omega_{\nu}}}-\log{\Pi_{\omega_{\nu-1}}},$ (25a)
$\displaystyle\sigma_{\mathrm{vis}}[\underline{\omega},\nu]$
$\displaystyle=-\log{p_{\omega_{\nu}}^{(\nu)}}+\log{p_{\omega_{\nu-1}}^{(\nu-1)}},$
(25b) $\displaystyle\sigma_{\mathrm{mot}}[\underline{\omega},\nu]$
$\displaystyle=B^{\omega_{\nu-1}}_{\omega_{\nu}}$ (25c)

can be used to construct the previously known fundamental forms of entropy
variation Schnakenberg (1976):

$\displaystyle\sigma_{\mathrm{sys}}[\underline{\omega},\nu]$
$\displaystyle=\sigma_{\mathrm{int}}[\underline{\omega},\nu]+\sigma_{\mathrm{vis}}[\underline{\omega},\nu],$
(26a) $\displaystyle\sigma_{\mathrm{tot}}[\underline{\omega},\nu]$
$\displaystyle=\sigma_{\mathrm{vis}}[\underline{\omega},\nu]+\sigma_{\mathrm{mot}}[\underline{\omega},\nu],$
(26b) $\displaystyle\sigma_{\mathrm{med}}[\underline{\omega},\nu]$
$\displaystyle=\sigma_{\mathrm{mot}}[\underline{\omega},\nu]-\sigma_{\mathrm{int}}[\underline{\omega},\nu].$
(26c)

The system entropy variation $\sigma_{\mathrm{sys}}$ is equivalent to the
change of coarse-grained entropy. It consists of the visible entropy of the
ensemble and the intrinsic entropy of the current state. Often, the latter is
forgotten Seifert (2011).

We can write the average of the _total entropy production_
$\sigma_{\mathrm{tot}}$ as a Kullback-Leiber divergence,

$\displaystyle\left\langle\\!\left\langle\sigma_{\mathrm{tot}}\right\rangle\\!\right\rangle^{(\nu)}=$
$\displaystyle\sum_{i,j}\left[p_{i}^{(\nu)}s^{i}_{j}\log\frac{p_{i}^{(\nu)}s^{i}_{j}}{p_{j}^{(\nu+1)}s^{j}_{i}}\right]\geq
0,$

which is always positive. It is the temporal variation of the difference
$\left\langle\\!\left\langle
s_{\mathrm{cg}}-s_{\mathrm{fg}}\right\rangle\\!\right\rangle$, which itself is
a Kullback-Leibler divergence. It characterizes the information gain of the
microscopic over the coarse-grained density. Its positive variation (i.e. its
monotonous growth) reflects the fact that one cannot resolve the ever-refining
features of the microscopic density.

The _entropy change in the medium_ $\sigma_{\mathrm{med}}$ consists of a
motance and an intrinsic part. It has been understood that the latter is
important for a consistent thermodynamic interpretation Seifert (2011, 2012).
From Eqs. (19) and (12) we see that it is directly related to phase-space
contraction:

$\sigma_{\mathrm{med}}[\underline{\omega}]=\log\frac{s^{\omega_{\nu-1}}_{\omega_{\nu}}}{s^{\omega_{\nu}}_{\omega_{\nu-1}}}-\log\frac{\Pi_{\omega_{\nu}}}{\Pi_{\omega_{\nu-1}}}=\log\eta^{\omega_{\nu-1}}_{\omega_{\nu}}$
(27)

For maps, the phase-space contraction per unit time, $\Lambda=\log J$ Ruelle
(1999); Rondoni and Mejia-Monasterio (2007), is given by the logarithm of the
Jacobian determinant
$J=\left|\frac{\,\mathrm{d}\Phi\,}{\,\mathrm{d}x\,}\right|$ of the dynamics.
In the current affine-linear case, $J|_{x}=\eta^{\omega_{k-1}}_{\omega_{k}}$
for $x\in\mathcal{C}^{\omega_{k-1}}_{\omega_{k}}$.

## Connection to the continuous description

For us, NMBM serve as a model for a _stroboscopic observation_ of physical
dynamics with a time lag $\tau_{\mathrm{obs}}$. However, stochastic
thermodynamics is usually formulated using _continuous-time Markov processes_
rather than Markov chains. In the continuous time limit, the Master equation
(15) becomes a differential equation
$\dot{p}_{i}(t)=\sum_{j}W^{j}_{i}p_{j}(t)$. Its infinitesimal generator
$W^{j}_{i}$ consists of transition rates $w^{i}_{j}$, which can be thought of
as the $\tau_{\mathrm{obs}}\to 0$ limit of transition probability divided by
$\tau_{\mathrm{obs}}$. The continuous-time versions of the fundamental entropy
and variation functionals have been identified previously Seifert (2005). They
resemble the discrete-time functionals with $s^{i}_{j}$ substituted by
$w^{i}_{j}$ and finite-time variations substituted by derivatives. However,
there are some subtleties to be taken care of. Firstly, in the continuous-time
description, the measure $\mathbb{P}[\underline{\omega}]$ of a finite
trajectory is zero and summing over trajectories involves integration over
stochastic jump times. Secondly, in the present case, the average value
$\left\langle\\!\left\langle\sigma\right\rangle\\!\right\rangle^{(\nu)}\equiv\left\langle\\!\left\langle
s[\underline{\omega},\nu]-s[\underline{\omega}^{(\nu-1)},\nu-1]\right\rangle\\!\right\rangle^{(\nu)}$
of a variation is the same as the change in one time-step of the averaged
entropies, $\left\langle\\!\left\langle
s\right\rangle\\!\right\rangle^{(\nu)}-\left\langle\\!\left\langle
s\right\rangle\\!\right\rangle^{\nu-1}$. This does not apply to the
continuous-time case, because averaging and taking the time-derivative do not
commute for non-stationary probabilities. Hence, temporal boundary-terms need
to be attributed for. Thirdly, we can rewrite the average total entropy
production in the following form:

$\displaystyle\left\langle\\!\left\langle\sigma_{\mathrm{tot}}\right\rangle\\!\right\rangle^{(\nu)}=$
$\displaystyle\sum_{i,j}\left[p_{i}^{(\nu)}s^{i}_{j}\log\frac{p_{i}^{(\nu)}s^{i}_{j}}{p_{j}^{(\nu)}s^{j}_{i}}\right]+\sum_{i}\left[p_{i}^{(\nu+1)}\log\frac{p_{i}^{(\nu)}}{p_{i}^{(\nu+1)}}\right]$

Substituting $s^{i}_{j}$ by $w^{i}_{j}$ in the first term, one obtains the
familiar term used in the continuous-time case, which is already a Kullback-
Leibler divergence on its own. However, for the discrete-time case the second
term is needed to account for the change of the ensemble between two jumps. It
vanishes for $\tau_{\mathrm{obs}}\to 0$ because $p_{i}^{(\nu)}\to
p_{i}^{(\nu+1)}$.

## Phase-space preserving systems

Finally, we demonstrate the consistency of our approach in the context of
phase-space preserving maps. From the _dynamical systems perspective_ , we
consider such maps as stroboscopic pictures of a Hamiltonian dynamics.
_Thermodynamically_ , a Hamiltonian system is an isolated system which does
not exchange energy, particles or entropy with its environment at any time.
Its steady-state distribution is _microcanonical_ , i.e.  uniform on phase
space. In _stochastic terms_ , the steady state of an isolated system features
_detailed balance_ Schnakenberg (1976) for the steady-state distribution
$p^{\infty}_{i}$:

$p_{i}^{\infty}s^{i}_{j}=p_{j}^{\infty}s^{j}_{i}.$ (28)

Henceforth, we assume Eq. (28) and reversibility, Eq. (4), to show that the
interpretations are consistent. Definition (14) directly yields
$s^{i}_{j}=c_{ij}\Pi_{j}$ with some constant $c_{ij}=c_{ji}$. Because
$\sum_{j}s^{i}_{j}=1,\,\forall i$, it follows that $c_{ji}=1$, $\forall i,j$.
With $p_{i}^{\prime}$ denoting the iterated values of the probabilities
$p_{i}$ in Eq. (15), we hence find

$\displaystyle\left\langle\\!\left\langle\sigma_{\mathrm{med}}\right\rangle\\!\right\rangle$
$\displaystyle=\sum_{i,j}\left[p_{i}s^{i}_{j}\log\frac{s^{i}_{j}}{s^{j}_{i}}\right]-\sum_{i}\left[\left(p^{\prime}_{i}-p_{i}\right)\log\Pi_{i}\right]$
$\displaystyle=\sum_{i,j}\left[p_{i}\Pi_{j}\log\Pi_{j}-p_{i}\Pi_{j}\log\Pi_{i}-p_{j}\Pi_{i}\log\Pi_{i}\right]$
$\displaystyle\quad+\sum_{i}\left[p_{i}\log\Pi_{i}\right]=\sum_{i}\left[\left(-p_{i}+p_{i})\log\Pi_{i}\right)\right]=0.$

Hence, the average entropy change in the medium—which equals the average
phase-space contraction—identically vanishes _for all times and independent of
the initial ensemble_ , as is expected.

## Summary and outlook

We have shown how for NMBM the functionals of stochastic thermodynamics emerge
from information theory without relying on thermodynamic assumptions. We
presented interpretations of the different notions of entropy and entropy
variation, which are consistent both with information theory and work on non-
phase-space preserving “thermostatted” dynamics.

In analogy to the role of Smale’s horseshoe Smale (1967) as a hallmark of
chaos, NMBM could hence help understanding salient features of reversible
chaotic dynamics and their connection to thermodynamics. A possible approach
in that direction is to treat the NMBM vertices as linearizations of
fundamental structures in such dynamical systems. Then, a coarse-graining
procedure under the assumption of separation of scales, in the spirit of Ref.
Esposito (2012) could yield a better dynamical picture of local equilibrium.

###### Acknowledgements.

The authors are indebted to A. Wachtel, J. Horowitz and L. Rondoni for
illuminating discussions.

## References

  * Seifert (2008) U. Seifert, Eur. Phys. J. B 64, 423 (2008).
  * Maes (2004) C. Maes, in  _Poincaré Seminar 2003: Bose-Einstein condensation-entropy_ (Birkhäuser, Basel, 2004) p. 145.
  * Seifert (2005) U. Seifert, Phys. Rev. Lett. 95, 40602 (2005).
  * Seifert (2012) U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
  * Kurchan (1998) J. Kurchan, J. Phys. A 31, 3719 (1998).
  * Lebowitz and Spohn (1999) J. Lebowitz and H. Spohn, J. Stat. Phys. 95, 333 (1999).
  * Jarzynski (1997) C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
  * Crooks (1999) G. Crooks, Phys. Rev. E. 60, 2721 (1999).
  * Esposito  _et al._ (2009) M. Esposito, K. Lindenberg, and C. Van den Broeck, Phys. Rev. Lett. 102, 130602 (2009).
  * Horowitz and Parrondo (2011) J. Horowitz and J. Parrondo, New. J. Phys. 13, 123019 (2011).
  * Sagawa and Ueda (2012) T. Sagawa and M. Ueda, Phys. Rev. E. 85, 021104 (2012).
  * Toyabe  _et al._ (2010) S. Toyabe, T. Sagawa, M. Ueda, E. Muneyuki, and M. Sano, Nature Physics (2010).
  * Tietz  _et al._ (2006) C. Tietz, S. Schuler, T. Speck, U. Seifert, and J. Wrachtrup, Phys. Rev. Lett. 97, 50602 (2006).
  * Ciliberto  _et al._ (2010) S. Ciliberto, S. Joubaud, and A. Petrosyan, J. Stat. Mech. 2010, P12003 (2010).
  * Seifert (2011) U. Seifert, Eur. Phys. J. E 34, 1 (2011).
  * Esposito (2012) M. Esposito, Phys. Rev. E. 85, 041125 (2012).
  * Speck and Seifert (2007) T. Speck and U. Seifert, J. Stat. Mech. 2007, L09002 (2007).
  * Vollmer (2002) J. Vollmer, Phys. Rep. 372, 131 (2002).
  * Colangeli  _et al._ (2011) M. Colangeli, R. Klages, P. De Gregorio, and L. Rondoni, J. Stat. Mech. 2011, P04021 (2011).
  * Hopf (1948) E. Hopf,  _Ergodentheorie_ (Chelsea Publishing Company, 1948).
  * Gaspard (2005) P. Gaspard,  _Chaos, scattering and statistical mechanics_ , Vol. 9 (Cambridge University Press, 2005).
  * Vollmer  _et al._ (1997) J. Vollmer, T. Tél, and W. Breymann, Phys. Rev. Lett. 79, 2759 (1997).
  * Rondoni  _et al._ (2000) L. Rondoni, T. Tél, and J. Vollmer, Phys. Rev. E. 61, 4679 (2000).
  * Evans and Searles (2002) D. Evans and D. Searles, Adv. Phys. 51, 1529 (2002).
  * Maes and Netočný (2003) C. Maes and K. Netočný, J. Stat. Phys. 110, 269 (2003).
  * Tel and Vollmer (2000) T. Tel and J. Vollmer, in  _Hard Ball Systems and the Lorentz Gas_ (Springer Berlin, 2000) pp. 367–418.
  * Beck and Schögl (1995) C. Beck and F. Schögl,  _Thermodynamics of chaotic systems: an introduction_ , Vol. 4 (Cambridge University Press, 1995).
  * Shannon (1948) C. E. Shannon, Bell System Technical Journal 27, 379 (1948), many reprints available.
  * Jepps and Rondoni (2010) O. Jepps and L. Rondoni, J. Phys. A 43, 133001 (2010).
  * Altaner  _et al._ (2012) B. Altaner, S. Grosskinsky, S. Herminghaus, L. Katthän, M. Timme, and J. Vollmer, Phys. Rev. E. 85, 041133 (2012).
  * Schnakenberg (1976) J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976).
  * Ruelle (1999) D. Ruelle, J. Stat. Phys. 95, 393 (1999).
  * Rondoni and Mejia-Monasterio (2007) L. Rondoni and C. Mejia-Monasterio, Nonlinearity 20, R1 (2007).
  * Smale (1967) S. Smale, Bull. Amer. Math. Soc 73, 747 (1967).


Classical and quantum gravity from relativistic quantum mechanics

Walter Smilga

Geretsried, Germany

<EMAIL_ADDRESS>
## Abstract

It is common practice to describe elementary particles by irreducible unitary
representations of the Poincaré group. In the same way, multi-particle systems
can be described by irreducible unitary representations of the Poincaré group.
Representations of the Poincaré group are characterised by fixed eigenvalues
of two Casimir operators corresponding to a fixed mass and a fixed angular
momentum. In multi-particle systems (of massive spinless particles), fixing
these eigenvalues leads to correlations between the particles. In the quasi-
classical approximation of large quantum numbers, these correlations take on
the structure of a gravitational interaction described by the field equations
of conformal gravity. A theoretical value of the corresponding gravitational
constant is calculated. It agrees with the empirical value used in the field
equations of general relativity.

##  1 Introduction

Although the axioms of (non-relativistic) quantum mechanics have long since
found a generally accepted form, their validity for relativistically
interacting multi-particle systems has not yet been seriously considered. The
general view is still that such systems must be described by local quantum
field theories111Werner Heisenberg in his last lecture on elementary particle
physics (1965): ‘Today we know that a theory of elementary particles must be a
local quantum field theory.’. This is why we now rely on a model based on
interacting quantised fields, the Standard Model of particle physics, which
describes three of the four fundamental interactions. Before the development
of non-perturbative methods (e.g. lattice QCD), numerical results from this
model were obtained using perturbation algorithms. Unfortunately, these
algorithms lead to mathematical inconsistencies at higher orders. These
inconsistencies have prevented a full understanding of the mathematical and
physical content of the model. Instead, they have opened the door to fanciful,
reifying interpretations [1] (‘virtual particles constantly popping in and out
of existence’222Scientific American, 9 October 2006), which have further
blocked the way to a rational understanding of the model. Whether the purely
numerical lattice methods will lead to a deeper understanding remains to be
seen.

If by understanding we mean a mathematically consistent traceability to the
two fundamental principles of nature, as formulated by the axioms of quantum
mechanics and Poincaré invariance, then we certainly have not yet understood
the Standard Model. After decades of stagnation on this issue, it seems that
only an unbiased, model-independent view of relativistic particle physics can
lead us out of this impasse.

A promising perspective emerges from the following observation: Just like
linking two angular momenta to a resulting angular momentum, linking two
single-particle states, each belonging to an irreducible representation of the
Poincaré group, to a physically meaningful two-particle state is by no means a
trivial task. It is not sufficient to simply form the product state of two
single-particle states; the representation theory of the Poincaré group rather
tells us, that the product representation obtained this way must be reduced to
a direct sum of irreducible representations in order to obtain descriptions of
two-particle systems with well-defined linear and angular momenta. The
resulting two-particle states are then obtained as superpositions of product
states with coefficients equal to the Clebsch–Gordan coefficients of the
Poincaré group [2]. The Standard Model, although based on product
representations, clearly ignores this rule.

This article examines the consequences of strict adherence to this rule.

##  2 Axioms of relativistic quantum mechanics

Dirac’s axioms of quantum mechanics [3] can be formulated as follows333from
lecture notes http://web.mit.edu/8.05/handouts/jaffe1.pdf:

First Postulate: At each instant the state of a physical system is represented
by a normalised ket $\left.\mbox{\textbar}\psi\right>$ in the space of states.

Second Postulate: Every observable attribute of a physical system is described
by an operator that acts on the kets that describe the system.

Third Postulate: The only possible result of the measurement of an observable
$\mathcal{A}$ is one of the eigenvalues of the corresponding operator
$\hat{\mathcal{A}}$.

Fourth Postulate: When a measurement of an observable $\mathcal{A}$ is made on
a given normalised state $\left.\mbox{\textbar}\psi\right>$, the probability
of obtaining an eigenvalue $a_{n}$ is given by the square of the inner product
of $\left.\mbox{\textbar}\psi\right>$ with the normalised eigenstate
$\left.\mbox{\textbar}a_{n}\right>$, $|{\left<a_{n}|\psi\right>}|^{2}$.

Fifth Postulate: Immediately after the measurement of an observable
$\mathcal{A}$ has yielded a value $a_{n}$, the state of the system is the
normalised eigenstate $\left.\mbox{\textbar}a_{n}\right>$.

Sixth Postulate: The time evolution of the state of a quantum system is
described by
$\left.\mbox{\textbar}\psi(t)\right>=\hat{U}(t,t_{0})\left.\mbox{\textbar}\psi(t_{0})\right>$,
for some unitary operator $\hat{U}$.

Note that the sixth axiom is only valid as long as no measurement is made
according to the fifth axiom. Ignoring this restriction leads leads to the so-
called ‘measurement problem’.

In this general form, these axioms do not distinguish between non-relativistic
and relativistic particle systems.

If $\left.\mbox{\textbar}\psi\right>$ belongs to a product representation of
the Poincaré group, it is a trivial exercise to write down the corresponding
operator $\hat{U}(t,t_{0})$ (see, e.g.,[4], pp. 148–150). Of course, a product
representation only describes ‘free’, i.e. non-interacting particles.

Historically, the challenge has been to add the four fundamental interactions
in a way that preserves relativistic covariance. This led to the Standard
Model of particle physics. Instead of relying on the results of the
representation theory of the Poincaré group [5], this model explicitly models
interactions by additional interaction-mediating fields. It redefines
$\left.\mbox{\textbar}\psi\right>$ as an operator-valued field based on a
Fock-space representation [6], misleadingly called ‘second quantisation’. When
applied to a ‘vacuum’ state $\left.\mbox{\textbar}0\right>$, this field
operator generates a multi-particle state. The operator $\hat{U}(t,t_{0})$ is
then derived from a Lagrangian based on the matter field (fermion field) and a
bosonic field that provides the interaction, both linked by an interaction
term. Transition amplitudes are determined using a perturbation algorithm in
which the interaction term is treated as a (small) perturbation of the matter
field.

Unfortunately, this perturbation algorithm leads to mathematical
inconsistencies that can only be solved by a formal trick called
‘renormalisation’. Since this trick does not work for the gravitational
interaction, the Standard Model describes only the electromagnetic, weak, and
strong interactions.

This article takes an alternative, non-perturbative approach. It is firmly
based on the representation theory of the Poincaré group. After Wigner’s
seminal paper in 1939 [5], elementary particles have been successfully
described by irreducible unitary representations of the Poincaré group. It is
therefore only logical to also describe multi-particle configurations of
elementary particles by irreducible unitary representations of the Poincaré
group and to consider $\left.\mbox{\textbar}\psi\right>$ as a state of such a
representation. Under this premise, the above axioms can be regarded as axioms
of relativistic quantum mechanics.

##  3 Irreducible representations of the Poincaré group

Irreducible unitary representations of the Poincaré group $P(3,1)$ are
characterised by fixed eigenvalues of two Casimir operators, corresponding to
a mass and an intrinsic angular momentum. This is analogous to the mass and
spin of irreducible single-particle representations (see textbooks on quantum
field theory, e.g.[4], pp. 36–53). The state space of such a representations
has an orthonormal basis consisting of eigenstates of a complete set of
commuting operators, which, in addition to the two Casimir operators, can be
the four components of the (total) momentum and a component of the
pseudovector $\frac{1}{2}\epsilon_{\sigma\mu\nu\lambda}M^{\mu\nu}p^{\lambda}$,
which in the rest frame corresponds to a component of the angular momentum
$M^{ik}$. In contrast to single-particle representations, where the spin
determines the angular momentum, in multi-particle representations an
intrinsic orbital angular momentum also contributes to the total angular
momentum.

The reduction of a product representation to a sum of irreducible
representations implies a decomposition of the associated state space, which
in turn leads to correlations between the individual eigenstates of momentum
$\left.\mbox{\textbar}\mathbf{p}_{1}\right>$ and
$\left.\mbox{\textbar}\mathbf{p}_{2}\right>$: Fixing the eigenvalue of the
first Casimir operator correlates the particle states in such a way that the
mass, defined by the total momentum
$\mathbf{p}=\mathbf{p}_{1}+\mathbf{p}_{2}$, takes on a fixed value. Fixing the
eigenvalue of the second Casimir operator forces the particle states into a
rotationally symmetric momentum-entangled structure, such as in the two-
particle eigenstate of the total linear and angular momentum (see [7])

$\left.\mbox{\textbar}\mathbf{p},w\right>=\int_{\Omega}\\!\mathrm{d}\omega_{\mathbf{k}}\left.\mbox{\textbar}\mathbf{p}_{1}+\mathbf{k}\right>c_{\mathbf{p},w}(\mathbf{k})\left.\mbox{\textbar}\mathbf{p}_{2}-\mathbf{k}\right>.$
(1)

The integration over $\mathbf{k}$ is over a circular path $\Omega$. This
structure can be understood as description of two particles correlated
(entangled) by the virtual exchange of momenta $\mathbf{k}$.

Integrals with a similar structure can be found in the perturbation algorithm
of the Standard Model. There, however, the integration is performed over the
entire three-dimensional momentum space. For loop-like Feynman diagrams, this
leads directly to the well-known divergences.

The two-particle state (1) clearly shows that the particle momenta cannot (not
even approximately) be eigenstates simultaneously with the total linear
momentum and angular momentum. This property already follows from the
commutation relations of the Poincaré group: In fact, the commutation
relations of the infinitesimal generators of $P(3,1)$ (see [4], p. 45) allow
the simultaneous diagonalisation of a linear momentum (the total momentum) and
of an angular momentum component with an axis of rotation parallel to the
linear momentum; they do not, however, allow the simultaneous diagonalisation
of another momentum that is not parallel to the axis of rotation.

The rotational symmetry of the eigenstates of the orbital angular momentum is
preserved in the quasi-classical limit of large quantum numbers. The
corresponding wave functions define probability distributions in the form of
well-defined rotationally symmetric closed orbits (see [8], pp. 27–29). Wave
packets travelling along such orbits necessarily have time-dependent momenta
that are subject to the conservation of the total momentum.

##  4 Conservation of total momentum in multi-particle systems

In the quasi-classical limit, the conserved total momentum of the multi-
particle system under consideration can therefore be written as a sum over
time-dependent particle momenta. In the rest frame of the system, the
conservation law of momentum takes the form

$\sum_{i=1}^{N}\mathbf{p}_{i}(t)=\mathbf{0}.$ (2)

The parameter $t$ denotes the time of a time cut.

In Newtonian mechanics, a time-dependent particle momentum $\mathbf{p}_{i}(t)$
defines a force that is equal to the time derivative of that momentum,
whatever the reason for the time dependence. Since the total momentum is
conserved, this force must be compensated by forces from other particles. Thus
we can describe the time dependence of the particle momenta in terms of an
exchange of forces or, equivalently, by a virtual (instantaneous) exchange of
momentum between the particles.

Note that ‘virtual exchange of momentum’ is used here as a convenient
descriptive phrase for describing the structure of Eq. (2) and its
consequences. It is not meant in the sense of a real exchange of momentum,
which would be subject to relativistic causality.

This description is valid for arbitrary particle momenta and masses. It is
therefore invariant under rotations and boost operations applied to the
individual particle momenta. Furthermore, it is invariant under conformal
scaling of the individual particle momenta, i.e. individual scaling of the
particle masses. These properties define the intrinsic symmetries of Eq. (2).

If Eq. (2) is rewritten in the form

$\mathbf{p}_{1}(t)=-\sum_{i=2}^{N}\mathbf{p}_{i}(t),$ (3)

it is then clear that the trajectory of particle $1$ is completely determined
by the time-dependent momenta of particles $2$ to $N$, together with the
position and momentum of particle $1$ at $t=0$. Had we not known that we were
dealing with the conservation of total momentum, we might have suspected that
the mere presence of neighbouring matter would bend the trajectory of particle
$1$. Alternatively, we could have argued that the presence of matter bends
spacetime along the trajectory of particle $1$. This clearly evokes
associations with the phenomenon of gravitation and suggests that particles
$2$ to $N$ generate a virtual gravitational field, expressed by the metric of
spacetime, which bends the trajectory of particle $1$ in such a way that the
total momentum of the system is the same for all times $t$.

Despite the virtual nature of this field, its structure is in principle
measurable by observing the trajectories of particle $1$. A structure that can
be measured should also be theoretically describable. Of course, this
description must be compatible with Einstein’s principle of relativity.

##  5 Relativistic field equations

In establishing the field equations of general relativity, Einstein required
that ‘the general laws of nature must be expressed by equations that are valid
for all coordinate systems, that is, that are covariant (generally covariant)
with respect to arbitrary substitutions’ ([9], p. 776, translation by the
author). We adopt this requirement and also the coordinate-independent
description of curved trajectories in four-dimensional spacetime with the
tools of differential geometry.

The concept of virtual momentum exchange is easily brought into a four-
dimensional form by replacing $\mathbf{p}_{i}(t)$ by a momentum density
$\mathbf{p}_{i}(x)$ and combining $\mathbf{p}_{i}(x)$ and the corresponding
energy density $p^{0}_{i}(x)$ into the four-vector $p_{i}(x)$. As required by
Einstein, $p_{i}(x)$ then becomes part of a covariant energy-momentum tensor
$T^{\mu\nu}(x)$ (of massive point-like particles, see, e.g.,[10]).

The variations of the momentum along a trajectory ‘generate’ virtual exchange
momenta. Since these exchange momenta are derived from the time-dependent
particle momenta alone, they are independent of the trace $T$ of $T^{\mu\nu}$
corresponding to the (time-constant) particle masses. At the same time, other
particles ‘absorb’ the exchange momenta, causing their trajectories to change.
These changes can again be described in terms of time-dependent momenta or,
preferably, in terms of the curvatures of the trajectories, which are easier
observed by a distant observer. Again, only the traceless part of $T^{\mu\nu}$
is involved, namely

$T^{\mu\nu}-\textstyle{\frac{1}{4}}g^{\mu\nu}\,T.$ (4)

This is consitent with the fact that the masses of the particles do not enter
into the conservation law of momentum, Eq. (2).

The exchange process will now be formulated as a law of nature (in the sense
of Einstein’s requirement); that is, as a covariant relation between the
traceless energy-momentum tensor (emitter) and a geometric tensor $W^{\mu\nu}$
(absorber), derived from a curvature tensor. Together, these tensors are
expected to form the covariant equation

$W^{\mu\nu}=\kappa\,(T^{\mu\nu}-\textstyle{\frac{1}{4}}g^{\mu\nu}\,T),$ (5)

where $\kappa$ is a constant needed for dimensional reasons. As a counterpart
to the traceless energy-momentum tensor, $W^{\mu\nu}$ must also be traceless.

To determine the explicit form of $W^{\mu\nu}$, we follow Mannheim [11] and
use the variational principle of least action $\delta I=0$. This requires a
suitably chosen action $I=I_{W}+I_{M}$, the sum of a geometric action $I_{W}$
and a material action $I_{M}$, which is invariant under all relevant
coordinate transformations, including local conformal scaling. The functional
variation of $I$ with respect to the metric tensor $g^{\mu\nu}$ should then
give the expected covariant relation between the curvature tensor and the
energy-momentum tensor.

The only geometric action that is invariant under local conformal scaling is
the Weyl action (see [11])

$I_{W}=-\textstyle{\frac{1}{4}}\int\\!\mathrm{d}^{4}x\,(-g)^{1/2}\,C_{\lambda\mu\nu\xi}C^{\lambda\mu\nu\xi},$
(6)

where $C_{\lambda\mu\nu\xi}$ is the conformally invariant Weyl tensor, which
is the traceless part of the Riemann curvature tensor $R_{\lambda\mu\nu\xi}$.

For the material action,

$I_{M}=\kappa\int\\!\\!\mathrm{d}^{4}x\,(-g)^{1/2}\,g_{\mu\nu}\,(T^{\mu\nu}-\textstyle{\frac{1}{4}}g^{\mu\nu}\,T)$
(7)

can be chosen. The term under the integral is the trace of the traceless
energy-momentum tensor, which is zero in any coordinate system:
$g_{\mu\nu}\,(T^{\mu\nu}-\textstyle{\frac{1}{4}}g^{\mu\nu}\,T)\ =\\\
g_{\mu\nu}\,T^{\mu\nu}-T=T-T=0$. Thus $I_{M}$ is trivially invariant under
local conformal scaling.

The factor $(-g)^{1/2}$ in the integrals $I_{W}$ and $I_{M}$ can be dropped if
$x$ refers to a fixed Minkowskian coordinate system (see [9], p. 801).

As shown in [11], the functional derivative of $I_{W}+I_{M}$ with respect to
$g^{\mu\nu}$ yields the fourth-order partial differential equations

$W^{\mu\nu}=\kappa\,(T^{\mu\nu}-\textstyle{\frac{1}{4}}g^{\mu\nu}\,T),$ (8)

where $W^{\mu\nu}$, expressed in terms of the Riemann tensor and Ricci tensor
$R^{\mu\nu}$, has the form

$\displaystyle W^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}g^{\mu\nu}(R^{\alpha}_{\;\;\alpha})^{;\beta}_{\;\;\;;\beta}+R^{\mu\nu;\beta}_{\;\;\;\;\;\;\;;\beta}-R^{\mu\beta;\nu}_{\;\;\;\;\;\;\;;\beta}-R^{\nu\beta;\mu}_{\;\;\;\;\;\;\;;\beta}-2R^{\mu\beta}R^{\nu}_{\;\;\beta}\hskip
48.36958pt$ (9) $\displaystyle+$
$\displaystyle\frac{1}{2}g^{\mu\nu}R_{\alpha\beta}R^{\alpha\beta}-\frac{2}{3}g^{\mu\nu}(R^{\alpha}_{\;\;\alpha})^{;\beta}_{\;\;\;;\beta}+\frac{2}{3}(R^{\alpha}_{\;\;\alpha})^{;\mu;\nu}+\frac{2}{3}R^{\alpha}_{\;\;\alpha}R^{\mu\nu}$
$\displaystyle-$
$\displaystyle\frac{1}{6}g^{\mu\nu}(R^{\alpha}_{\;\;\alpha})^{2}$
$\displaystyle=$ $\displaystyle
2\,C^{\mu\lambda\nu\xi}_{\;\;\;\;\;\;\;\;\;\;\;;\lambda;\xi}-C^{\mu\lambda\nu\xi}R_{\lambda\xi}.$

Due to their derivation by a variational principle, these field equations
describe an equilibrium between emitted and absorbed virtual momenta.
Surprisingly, these equations have the same structure as Mannheim’s field
equations of conformal gravity, although the latter are based on completely
different physical premises. While Mannheim postulated a gauge invariance of
nature under local conformal scaling, conformal covariance here refers to the
aforementioned property of the conservation law of momentum to be covariant
under particle-individual conformal scaling. In contrast to Mannheim’s
approach, the particles involved here need not be massless.

As in Einstein’s general theory of relativity, the metric tensor
$g^{\mu\nu}(x)$ defines a pseudo-Riemannian manifold. Observers moving freely
along a trajectory in this manifold will experience a curved spacetime,
whereas they will infer the existence of a gravitational field from a
Minkowskian coordinate system.

Mathematically, the field equations are inhomogeneous differential equations
with the traceless energy-momentum tensor as the source (of virtual momenta).
They exploit the equivalence between the description of trajectories in terms
of time-dependent momenta – that is, by
$T^{\mu\nu}-\textstyle{\frac{1}{4}}g^{\mu\nu}\,T$ – and their description in
terms of curvatures – that is, by $W^{\mu\nu}$.

Since these are differential equations, a virtual momentum generated by a
particle at position $x$ can be absorbed by a second particle at or near
position $x$, but also by particles at other positions. Therefore, only a
fraction of the virtual momentum will be absorbed by the second particle.
Since Eq. (2) does not favour any particle, this must also apply to the field
equations. Quantum mechanically, this means that all particles have the same
initial probability of being involved in the exchange process. Therefore,
$\kappa$ must contain a constant probability factor of $1/N$, where $N$ is the
total number of particles. This gives $\kappa$ the role of a coupling constant
(see next section).

When applied to a real multi-particle system (subject to the conservation law
of momentum), the field equations translate the dynamical constraints on the
particle momenta, i.e. the conservation of total momentum expressed by Eq.
(2), into equivalent kinematic constraints on the particle trajectories in
spacetime. Thus, the field equations describe the structure of spacetime of
isolated multi-particle systems when the energy-momentum tensor is given, and
vice versa. As Wheeler put it in [12] (slightly modified by the author):
Matter tells spacetime how to curve; spacetime tells matter how to move.

##  6 Coupling constants

The constant $\kappa$ in the field equations (8) corresponds to the
gravitational constant in Einstein’s field equations. It connects the energy-
momentum tensor, which has the dimension of an energy density, i.e.
$\mathbf{J\,m^{-3}}$ or equivalently $\mathbf{kg\,m^{-1}s^{-2}}$, with a
squared curvature, which has the dimension $\mathbf{m^{-4}}$. Therefore,
$\kappa$ has the dimension $\mathbf{kg^{-1}m^{-3}s^{2}}$ or equivalently
$\mathbf{J^{-1}\,m^{-1}}$. It converts the dimension of the energy into the
dimension of an inverse length, which can be understood as the dimension of
the curvature of a trajectory.

Remarkably, there are two constants of nature, which together have the same
dimension as $\kappa$: The Planck constant $h$ with the dimension
$\mathbf{J\,s}$ and the speed of light $c$ with the dimension
$\mathbf{m\,s^{-1}}$ form the constant $1/(h\,c)$, which has the dimension
$\mathbf{J^{-1}\,m^{-1}}$. Its numerical value is $1/(6.626\times
10^{-34}\times 2.998\times 10^{8})$ $=5.034\times 10^{24}$.

Together with the dimensionless factor $1/N$, where $N$ is tentatively
$2.4\times 10^{67}$, the estimated number of atoms in our galaxy444Steve Cavil
on Oxford Education Board: “Our galaxy, the Milky Way, contains approximately
100 to 400 billion stars. If we take this as 200 billion or $2\times 10^{11}$
stars and assume that our sun is a reasonable average size we can calculate
that our galaxy contains about $(1.2\times 10^{56})\times(2\times
10^{11})=2.4\times 10^{67}$ atoms.”, we get an estimate of $\kappa$:

$\kappa=1/(N\,h\,c)\approx 2.1\times 10^{-43}\;\mathbf{kg^{-1}m^{-3}s^{2}}.$
(10)

For our solar system the equations of conformal gravity are equivalent to
those of general relativity (see [11]). This estimate can therefore be
compared to the value of the gravitational constant in Einstein’s field
equations. Based on the measured value of Newton’s gravitational constant
$G=6.674\times 10^{-11}\;\mathbf{kg^{-1}m^{3}s^{-2}}$, this value is

$\kappa_{E}=8\pi G/c^{4}\approx 2.077\times
10^{-43}\;\;\mathbf{kg^{-1}m^{-1}s^{2}}.$ (11)

The fact that $\kappa$ can be expressed in terms of $h$, $c$, and the
dimensionless constant $N$ with a unique origin removes the motivation to
speculate about physics at Planck scales. The factor $1/N$ makes clear that
the gravitational interaction, unlike the electromagnetic interaction, is a
multi-particle interaction, which has an extremely large range due to the lack
of shielding. This factor also gives a simple answer to the so-called
hierarchy problem.

For the electromagnetic interaction, the physically relevant representation is
an irreducible two-particle representation, since in larger systems positively
and negatively charged particles tend to combine into a configuration that is
electrically neutral to the outside.

In the article [7] cited above, the structure of irreducible (massive) two-
particle representations was analysed leading to the correct value for the
electromagnetic fine-structure constant $\alpha$ (= 1/137.03608245 vs. CODATA
value $1/137.035999139$). The agreement of the calculated value with the
experimental value identifies the kinematic constraints of irreducible two-
particle representations as ‘electromagnetic interaction’ and thus provides
further support for the concept of irreducible multi-particle representations
of the Poincaré group.

##  7 Quantisation of spacetime

As described above, irreducible multi-particle representations of the Poincaré
group in the limit of large quantum numbers lead to curved trajectories in a
pseudo-Riemannian spacetime. The inversion of this limit can therefore be
understood as ‘quantisation of spacetime’. This means that continuous
spacetime is replaced by discrete trajectories, which are then replaced by
wave functions of an irreducible multi-particle representation of $P(3,1)$.
This quantisation rule is transparent and well founded, both physically and
mathematically. Unlike other attempts to quantise spacetime, such as loop
quantum gravity or string theory, it provides a quantisation of spacetime at
experimentally accessible scales.

##  8 Conclusions

This article provides a consistent foundation for the phenomenon of gravity in
both the quantum mechanical and classical domains. The foundation is based
solely on the first principles of quantum mechanics and Poincaré invariance,
unified in the representation theory of the Poincaré group.

In multi-particle systems, which are described quantum mechanically by
irreducible multi-particle representations of the Poincaré group, the
kinematics of the individual particles is constrained in the sense that the
individual particle states are superposed in such a way that they form a
common rotationally symmetric eigenstate of the orbital angular momentum. As a
consequence, in the quasi-classical limit, the trajectories of the particles
are not rectilinear, but form a pseudo-Riemannian manifold, which is described
by the field equations of conformal gravity (CG).

The field equations of CG are currently being discussed as an alternative to
the field equations of general relativity (GR) (see e.g. [11]). In the form
presented here, they apply to isolated multi-particle systems consisting of
massive spinless particles or macroscopic bodies. For the solar system, they
not only reproduce the Schwarzschild solutions of GR but (in contrast to GR)
they are also able to correctly describe galactic rotation curves without the
aid of exotic dark matter (see [11] and the references given there). Thus, not
only Einstein–Newton gravity but also the puzzling galactic rotation curves
can be regarded as manifestations of the constrained kinematics of isolated
multi-particle systems.

Applied to our galaxy, the derived field equations provide a gravitational
constant of the right order of magnitude. The numerical agreement of this
constant with Einstein’s gravitational constant confirms that the strength of
gravity in our solar system is indeed determined by the surrounding Milky Way
galaxy. This can be generalised to the statement that the gravitational
constant within a galaxy is specific to that galaxy.

If we agree to define quantum gravity as the description of multi-particle
systems that is consistent with the axioms of quantum mechanics and leads in
the quasi-classical limit to the field equations of CG, then irreducible
unitary multi-particle representations of the Poincaré group are the answer to
the current search for a quantum theory of gravity. These representations are
in principle well-understood and provide the simplest description that fits
this definition.

In summary, the above considerations show that it is possible to describe
relativistically interacting multi-particle systems on the basis of the axioms
of quantum mechanics, without having to resort to the problematic methods of
interacting quantum fields. The absence of explicit interaction terms means
that the mathematical complexity is drastically reduced compared to
perturbation schemes based on quantum field theories as known from the
Standard Model. A reformulation of the Standard Model based on irreducible
multi-particle representations of the Poincaré group would not only allow a
smooth integration of gravity into the model, but would also eliminate
renormalisation problems. The main advance over the current Standard Model is
that this approach not only avoids the inconsistencies of local field
theories, but also provides exactly the two known long-range interactions,
together with their coupling constants.

Of course, these results do not devalue the Standard Model as a collection of
powerful algorithms ready for practical computation. But they do call into
question the role of gauge invariance as the cause of interactions. And they
expose the myth of ‘virtual particles constantly popping in and out of
existence’ as merely an inappropriately reifying interpretation of the Fock
space-based perturbation algorithm or, in the words of A. N. Whitehead, as a
‘fallacy of misplaced concreteness’ [13].

### Funding information

No funds, grants, or other support were received.

## References

  * [1] D. Mermin, What’s bad about this habit?, Physics Today 62(5), 8 (2009); doi:10.1063/1.3141952
  * [2] H. Joos, Zur Darstellungstheorie der inhomogenen Lorentzgruppe als Grundlage quantenmechanischer Kinematik, Fortschritte der Physik 10(2), 65–146 (1962); doi:10.1002/prop.2180100302.
  * [3] P. A. M. Dirac, The Principles of Quantum Mechanics, 4th edn, Oxford University Press, Oxford (1958).
  * [4] S. S. Schweber, An Introduction to Relativistic Quantum Field Theory, 2nd edn, Harper & Row, New York (1962).
  * [5] E. Wigner, On Unitary Representations of the Inhomogeneous Lorentz Group, Ann. Math. 40(1), 149 (1939); doi:10.2307/1968551.
  * [6] V. Fock, Konfigurationsraum und zweite Quantelung, Z. Physik 75, 622 (1932); doi:10.1007/bf01344458.
  * [7] W. Smilga, Structural properties of irreducible two-particle representations of the Poincaré group, Found. Phys. 49(7), 728 (2019); doi:10.1007/s10701-019-00277-9.
  * [8] A. R. Edmonds, Angular Momentum in Quantum Mechanics (Investigations in Physics, vol 4), Princeton University Press (1957).
  * [9] A. Einstein, Die Grundlage der allgemeinen Relativitätstheorie, Ann. Phys. 49(7), 769 (1916); doi:10.1002/andp.19163540702.
  * [10] G. Ares de Parga and R. E. Gonzalez-Narvaez, A further analysis of the energy-momentum for a system of point particles and for a perfect fluid in finite volume, J. Phys.: Conf. Ser. 845, 012002 (2017); doi:10.1088/1742-6596/845/1/012002.
  * [11] P. D. Mannheim, Alternatives to dark matter and dark energy, Progr. Part. Nucl. Phys. 56(2), 340 (2006); doi:10.1016/j.ppnp.2005.08.001, arXiv:astro-ph/0505266v2.pdf.
  * [12] J. A. Wheeler and K. W. Ford, Geons, Black Holes and Quantum Foam: A Life in Physics, Norton, New York (1998).
  * [13] A. N. Whitehead, Science and the Modern World. Lowell Lectures, 1925, Pelican Mentor Books, New York (1948) (First edition: 1925).