# A dynamic $p$-Laplacian

Alvaro de Diego Unanue  Department of Mathematics, TUM School of Computation,
Information and Technology, Technical University of Munich, 80333 Munich,
Germany  Gary Froyland  School of Mathematics and Statistics, University of
New South Wales, Sydney NSW 2052, Australia  Oliver Junge  Department of
Mathematics, TUM School of Computation, Information and Technology, Technical
University of Munich, 80333 Munich, Germany  Péter Koltai  Department of
Mathematics, University of Bayreuth, 95440 Bayreuth, Germany

###### Abstract

We generalise the dynamic Laplacian introduced in (Froyland, 2015) to a
dynamic $p$-Laplacian, in analogy to the generalisation of the standard
$2$-Laplacian to the standard $p$-Laplacian for $p>1$. Spectral properties of
the dynamic Laplacian are connected to the geometric problem of finding
“coherent” sets with persistently small boundaries under dynamical evolution,
and we show that the dynamic $p$-Laplacian shares similar geometric
connections. In particular, we prove that the first eigenvalue of the dynamic
$p$-Laplacian with Dirichlet boundary conditions exists and converges to a
dynamic version of the Cheeger constant introduced in (Froyland, 2015) as
$p\rightarrow 1$. We develop a numerical scheme to estimate the leading
eigenfunctions of the (nonlinear) dynamic $p$-Laplacian, and through a series
of examples we investigate the behaviour of the level sets of these
eigenfunctions. These level sets define the boundaries of sets in the domain
of the dynamics that remain coherent under the dynamical evolution.

##  1 Introduction

Time-dependent or non-autonomous dynamical systems frequently exhibit complex,
turbulent motion; prominent examples are geophysical flows in the atmosphere
and ocean, industrial fluid mixers, and the dynamics of biofluids. These flows
often display heterogeneous behaviour, i.e. there exist regions of fluid that
are being _coherently_ transported and slowly dispersed relative to the
surrounding fluid. There is a large literature on the identification of
coherent structures of one sort or another. Our interest is primarily in
_geometric_ characterisations of mixing.

In [14], Froyland proposed to characterise a _coherent set_ by a small average
_evolved_ boundary area to volume ratio. If this ratio is high for a
particular region of fluid, it is likely that this fluid region becomes highly
filamented over a finite flow time and therefore is strongly susceptible to
small-scale diffusive processes. Indeed, tracking the growth in irregularity
of fluid interfaces is well-established as a means of quantifying mixing
properties of fluids [3, 35, 18, 4]. Such a characterisation extends classical
(non-dynamic) notions from isoperimetric theory, whereby the _Cheeger ratio_
of a set is the ratio of its boundary area (boundary codimension-one volume)
to its enclosed volume; formally, the infimand on the right-hand side of (1).
In the characterisation of [14], coherent sets have a small _dynamic Cheeger
ratio_ , defined formally in (5). The problem of finding sets with small
dynamic Cheeger ratio can be transformed into finding eigenfunctions of a
certain Laplace–Beltrami operator, called _dynamic Laplacian_ [14].

To set up our contributions in this work, we recap some classical
isoperimetric theory. Throughout the paper (except for one of the examples),
we consider a compact $d$-dimensional manifold $M\subset\mathbb{R}^{d}$ with
nonempty Lipschitz boundary. Classically, the _Cheeger constant_ of $M$ is
defined as [9]

$h(M):=\inf_{A\subset M}\frac{\ell_{d-1}(\partial A)}{\ell_{d}(A)},$ (1)

where the infimum is taken over all open submanifolds $A$ of $M$ with
$C^{\infty}$ boundary and compact closure and $\ell_{d}$ denotes
$d$-dimensional volume.

If a set $A$ has a Cheeger ratio of $h(M)$, it is called a _Cheeger set_. Note
that a Cheeger set does not have to have $C^{\infty}$ boundary and might
intersect the boundary of $M$.

The Federer–Fleming theorem [12, 25] provides a tight connection between
$h(M)$ and the functional _Sobolev constant_ defined by an infimum taken over
$C^{\infty}$ functions $u:M\to\mathbb{R}$ with compact support111For a vector-
valued function $v$, we denote by $\|v\|_{p}$ the $L^{p}$ norm of the
pointwise Euclidean length $|v|$ of $v$.

$s(M):=\inf_{u\in C^{\infty}_{c}(M)}\frac{\|\nabla u\|_{1}}{\|u\|_{1}},$ (2)

namely $h(M)=s(M)$ [22, Remark 2.1]. One may think of the supports of an
infimising sequence of functions $u$ in (2) defining an infimising sequence of
sets $A$ in (1). Replacing the 1-norm by the 2-norm in (2) produces the
related minimisation problem

$\lambda_{2}=\inf_{u\in C^{\infty}_{c}(M)}\frac{\|\nabla
u\|_{2}^{2}}{\|u\|_{2}^{2}}.$ (3)

In fact, $\lambda_{2}>0$ is the leading eigenvalue222Later we will denote by
$\lambda_{p}^{(k)}$ the $k^{\rm th}$ eigenvalue of the $p$-Laplacian. For now,
we suppress the superscript when $k=1$. of the (negative) 2-Laplace operator
on $M$ with homogeneous Dirichlet boundary conditions, i.e. $u=0$ on $\partial
M$ [8, Rayleigh’s Theorem]. A connection between the two problems is given by
the well-known _Cheeger inequality_

$h(M)\leq 2\sqrt{\lambda_{2}}\,.$ (4)

Informally one may think of the Laplace eigenfunction $u_{2}$ corresponding to
$\lambda_{2}$ as a smoothed version of the limit of an infimising sequence of
functions $u$ in (2), and one of the superlevel sets of $u_{2}$ as a smoothed
version of the limit of an infimising sequence of sets $A$ in (1).

Dynamic counterparts of these constructions (for the case of Neumann boundary
conditions) have been developed in [14], and are recalled in section 2. The
particular expressions for the Dirichlet boundary condition case, (6)–(9), may
be found in [16]. As in the static case, a suitable superlevel set of the
leading eigenfunction of the dynamic Laplacian yields a smoothed version of
the set of interest, namely a set with small averaged isoperimetric ratio,
i.e. a _coherent set_. Multiple coherent sets can be found directly from
suitable superlevel sets of eigenfunctions at higher eigenvalues of the
dynamic Laplacian [14, 7, 16, 17]. Computationally, the eigenproblem for the
dynamic Laplacian can efficiently be solved by the finite element method [16],
leading to a scheme for which only trajectory data is required.

While the use of the eigenfunctions is computationally advantageous on the one
hand, it is not _a priori_ clear how well the level sets of these
eigenfunctions approximate solutions of the original problem (1). On the other
hand, while it is desirable to directly solve the original formulation using
the 1-norm, the solutions to this problem are in general not smooth and
typically certain regularisations of the problem have to be considered [13].
Instead of regularising the $p=1$ eigenproblem, one can try to solve a
(nonlinear) eigenproblem for a $p$-norm with $1\lesssim p\leq 2$; this is the
approach we pursue here. We show that the 2-norm constructions of [14] can be
generalised to $p$-norm versions, obtaining a $p$-norm version of the dynamic
Sobolev constant and a dynamic $p$-Laplacian. We prove a dynamic $p$-Cheeger
inequality for the leading eigenvalue of the (nonlinear) dynamic $p$-Laplacian
(Theorem 22), generalising the dynamic Cheeger inequality of [14]. In Theorem
3.5 we show that in the limit $p\to 1$ this dynamic $p$-Cheeger inequality
interpolates to the dynamic Federer–Fleming equality of [14]. Thus, for $p$
close to 1, one expects most level sets of the dynamic $p$-Laplacian to have a
small dynamic Cheeger ratio (relative to other sets) and therefore to
represent a relatively coherent set.

The numerical computation of the eigenfunctions of the (dynamic) $p$-Laplacian
becomes more challenging as $p$ comes closer to $1$. We employ a scheme
proposed by Yao and Zhou [36] for a certain general class of nonlinear
eigenproblems. Our numerical experiments in section 5 allow us to compare
coherent sets that result from superlevel sets of eigenfunctions of the
dynamic $p$-Laplacian for $p$ close to 1 with those for $p=2$. In all cases,
we find the dynamic $p$-Laplacian eigenfunctions have profiles closer to
indicator functions than the smooth profile of the dynamic $2$-Laplacian.
Despite this, the difference between the level sets of the eigenfunctions that
minimise their Cheeger ratios in the $1\lesssim p$ and $p=2$ cases is small,
indicating that the coherent sets obtained through level sets of the dynamic
2-Laplacian are indeed accurate estimates of solutions of the dynamic Cheeger
problem.

##  2 Preliminaries

###  2.1 The dynamic isoperimetric problem

As alluded to in the introduction, dynamic versions of (1) and (2) can be
constructed [14] (see [16] for the case of Dirichlet boundary data as used
below). The underlying dynamical system is defined by a volume-preserving
diffeomorphism333In (9) below, from a classical point of view, we need to
assume $T$ to be a $C^{2}$-diffeomorphism. Since we always view $\Delta$ as
mapping into a dual space (see (10) and for the dynamic version (19)),
however, $C^{1}$ regularity of $T$ is actually sufficient. $T:M\rightarrow M$,
which in applications is often given by the flow map of a conservative vector
field over some finite time span. We denote by $T_{*}:u\mapsto u\circ T^{-1}$,
$u\in L^{p}(M)$, the _pushforward operator_ corresponding to $T$. Since we
assume $T$ to be volume-preserving, $T_{*}$ coincides with the _transfer_ (or
_Perron–Frobenius_) operator of $T$. Its adjoint is $T^{*}$, the Koopman
operator, with $T^{*}:u\mapsto u\circ T$, $u\in L^{q}(M)$.

For a single application of $T$, the _dynamic Cheeger ratio_ of a subset
$A\subset M$ is

$\frac{\ell_{d-1}(\partial A)+\ell_{d-1}(T(\partial A))}{2\ell_{d}(A)},$ (5)

and the _dynamic Cheeger constant_ of $M$ is

${h}^{D}(M,T):=\inf_{A\subset M}\frac{\ell_{d-1}(\partial
A)+\ell_{d-1}(T(\partial A))}{2\ell_{d}(A)},$ (6)

where the infimum is taken over open submanifolds $A$ with compact closure and
$C^{\infty}$ boundary. Similarly, the _dynamic Sobolev constant_ is

${s}^{D}(M,T):=\inf_{u\in C^{\infty}_{c}(M)}\frac{\|\nabla u\|_{1}+\|\nabla
T_{*}u\|_{1}}{2\|u\|_{1}}.$ (7)

As in the static case, one has a dynamic Federer–Fleming theorem, i.e.
equality of these two quantities ([14] for the Neumann case, [16] for the
Dirichlet case): ${h}^{D}(M,T)={s}^{D}(M,T)$. Modifying the norm in problem
(7) as in the static case yields the minimisation problem

${\lambda}^{D}_{2}=\inf_{u\in C^{\infty}_{c}(M)}\frac{\|\nabla
u\|_{2}^{2}+\|\nabla T_{*}u\|_{2}^{2}}{2\|u\|_{2}^{2}}.$ (8)

Also like in the static case, ${\lambda}^{D}_{2}>0$ is the leading eigenvalue
of the (negative) _dynamic Laplace operator_ [14]

${\Delta}^{D}:=\frac{1}{2}\left(\Delta+T^{*}\Delta T_{*}\right)$ (9)

on $M$ subject to homogeneous Dirichlet boundary conditions [16]. A suitable
superlevel set of an eigenfunction of ${\Delta}^{D}$ at ${\lambda}^{D}_{2}$
yields a smoothed version of the limit of an infimising sequence of sets $A$
in (6), i.e. a coherent set.

###  2.2 1-norm minimisation versus 2-norm minimisation

In order to illustrate the effect of the transition from the $1$-norm to the
$2$-norm in (2), we consider a minimiser of the static problem on the unit
square $M=[0,1]^{2}$ (more precisely, the limit of an infimising sequence). It
can be shown (see e.g. [20]) that almost all of its level sets coincide with
the minimiser of the Cheeger problem. For this simple domain, one can
explicitly calculate the Cheeger constant and the Cheeger set (e.g. [19]):

$h([0,1]^{2})=\frac{4-\pi}{2-\sqrt{\pi}}\approx 3.772.$

The Cheeger set is a square with rounded corners of radius
$1/h([0,1]^{2})\approx 0.265$ touching the boundary; see e.g. [20, Remark 5]
and Figure 1.

Figure 1: Left: minimiser of (2). Right: corresponding minimiser of (1). The
red curve delimits the Cheeger set.

In contrast, if we alter the $1$-norm to the $2$-norm, the corresponding
minimiser $u(x,y)=\sin(\pi xy)$ is quite smooth and its level sets vary
strongly. Nevertheless, at least in this simple example, the optimal
superlevel set – with boundary shown in red in Figure 2 – is close to the
minimiser of (1). Moreover, the Cheeger ratio of this superlevel set is
approximately $3.890$, which is close to the exact minimal value of around
$3.772$.

Figure 2: Graph of $u(x,y)=\sin(\pi xy)$, the minimiser of $\frac{\|\nabla
u\|_{2}^{2}}{\|u\|_{2}^{2}}$, $u|_{\partial M}=0$ (left) and the level sets
(right). The red curve is the level set of $u$ with optimal Cheeger ratio.

Our aim in the sequel is to consider the norms $\|\cdot\|_{p}$ as $p\in(1,2)$
and to investigate what happens for $p\rightarrow 1$. Ideally, the minimisers
become flatter away from the optimal level set, while the level sets get
closer to the optimum level set, so that the optimal level sets can be more
robustly extracted. We also anticipate that the infimum of the Rayleigh
quotient with respect to the $p$-norms converges to $h(M)$. In the classical
(static) setting the above behaviour of the level sets and the convergence of
the Rayleigh quotient can be proven (see e.g. [20]). Further, the “flatness”
of the eigenfunctions for $p$ close to $1$ can be observed in numerical
solutions (see e.g. [19, Figure 8]). However, no such attempt has been made so
far for a functional involving dynamics.

###  2.3 The $p$-Laplacian

In this section we recall results on the consequences of replacing the
$1$-norm by the $p$-norm in (2). The minimiser can then be found by computing
eigenfunctions of the so-called _$p$ -Laplacian_ (see e.g. the introduction of
[24]). In the limit $p\rightarrow 1$, these eigenfunctions allow us to recover
solutions of the Cheeger problem (1) (see 2.4 below). We will give a short
overview of the definition of the $p$-Laplacian and some known results that
link it to the Cheeger problem. In what follows we always assume that $1<p<2$
and denote by $q$ the conjugate of $p$ satisfying $\smash{1/p+1/q=1}$.

Let $W^{1,p}_{0}(M)$ denote the classical Sobolev space of once weakly
differentiable $L^{p}$ functions on $M$ that are zero on the boundary and have
weak first derivatives in $L^{p}$. Consider the functional
$F:W_{0}^{1,p}(M)\to\mathbb{R}$,

$F(u):=\|\nabla u\|_{p}^{p}.$

Its (Fréchet) derivative $F^{\prime}(u)$ at some $u\in W^{1,p}_{0}(M)$ is an
element of the dual space of $W^{1,p}_{0}(M)$, which is denoted by
$W^{1,p}_{0}(M)^{*}$. The action of $F^{\prime}(u)$ on a function $v\in
W^{1,p}_{0}(M)$ can be calculated by $F^{\prime}(u)v=\lim_{h\to
0}\frac{1}{h}(F(u+hv)-F(u))$ to be

$F^{\prime}(u)v=p\int_{M}\nabla u|\nabla u|^{p-2}\nabla v.$

The (classical) $p$-Laplacian $\Delta_{p}:W_{0}^{1,p}(M)\rightarrow
W^{1,p}_{0}(M)^{*}$ is defined as the nonlinear operator
${\Delta_{p}u:=-\frac{1}{p}F^{\prime}(u)}$, thus

$(\Delta_{p}u)v=-\int_{M}\nabla u|\nabla u|^{p-2}\nabla v.$ (10)

If $u$ is sufficiently regular, then the functional $\Delta_{p}u$ can even be
represented by a function; in a similar fashion to Riesz duality between
$L^{p}$ and $L^{q}$. Thus, identifying the integral operator with its integral
kernel, we sometimes write

$\Delta_{p}u=\mathrm{div}(|\nabla u|^{p-2}\nabla u)$

for $u$ regular enough (meaning the right-hand side is in $L^{q}(M)$). This
“strong form” of the operator $\Delta_{p}$ will be intuitive for the
eigenvalue problems we introduce below. For rigorous derivations we always
consider the “weak form” (10). For $p=2$, the $p$-Laplacian coincides with the
standard Laplacian.

We will mainly be interested in the variational problem

$\min F(u)\quad\text{such that}\quad G(u):=\|u\|_{p}^{p}=1.$ (11)

The Euler–Lagrange equation of (11) is

$F^{\prime}(u)=\lambda G^{\prime}(u).$

As the derivative of $G$ is $G^{\prime}(u)v=p\int_{M}u|u|^{p-2}v$, the
expression immediately above becomes

$\int_{M}\nabla u|\nabla u|^{p-2}\cdot\nabla
v=\lambda\int_{M}u|u|^{p-2}v\quad\text{for all }v\in W_{0}^{1,p}(M),$ (12)

or in its strong version

$-\Delta_{p}u=\lambda u|u|^{p-2}.$ (13)

If a pair $(\lambda,u)\in\mathbb{R}\times W_{0}^{1,p}(M)$ solves (12), we call
$\lambda$ an eigenvalue of $-\Delta_{p}$ and $u$ an eigenfunction. In contrast
to the linear case $p=2$, where a plethora of results is available, the set of
eigenvalues of the $p$-Laplacian is not fully known [28, 24]. One can
construct a subset of it by a min-max-principle based on the Krasnoselskii
genus [28], but it is not known whether these exhaust the spectrum of
$\Delta_{p}$.

In our case, we will be mostly concerned with the smallest eigenvalue, which
corresponds to the minimum of (11). Since $F$ and $G$ have the same order $p$
of homogeneity, one can equivalently search for a minimum of the functional
$J:W_{0}^{1,p}(M)\backslash\\{0\\}\to\mathbb{R}$,

$J(u)=\frac{F(u)}{G(u)}.$

######  Theorem 2.1.

[24, Lemma 5] Let $d\geq 2$, $1<p<2$ and $M\subseteq\mathbb{R}^{d}$ a bounded
manifold of dimension $d$ with Lipschitz boundary. Then

$\lambda_{p}:=\inf_{u\in W_{0}^{1,p}(M)\setminus\\{0\\}}J(u)$

is positive and achieved by a function $u_{p}\in W_{0}^{1,p}(M)$ satisfying
(12) with $\lambda=\lambda_{p}$.

######  Theorem 2.2.

[24, Lemma 5] Let $d\geq 2$, $1<p<2$ and $M\subset\mathbb{R}^{d}$ a bounded
manifold of dimension $d$ with Lipschitz boundary. If $u\in W^{1,p}_{0}(M)$ is
a solution of (11) then $u$ is also a solution of (12). Further, a solution
$u$ of (11) either fulfils $u>0$ or $u<0$ in the interior of $M$, i.e. it does
not change sign.

######  Theorem 2.3.

[23, Theorem 6] Let $d\geq 2$, $1<p<2$ and $M\subset\mathbb{R}^{d}$ a bounded
manifold of dimension $d$ with Lipschitz boundary and
$F(u):=\frac{1}{p}\|\nabla u\|_{p}^{p}$. Then all minima of $F$ are scalar
multiples of each other.

The preceding theorems show that there is a unique non-negative normalised
solution to (11), satisfying (12) with $\lambda=\lambda_{p}$. We call
$\lambda_{p}$ the _first eigenvalue_ of $-\Delta_{p}$ and $u_{p}$ the _first
eigenfunction_ of $-\Delta_{p}$ on $M$.

We are particularly interested in the following properties that link $u_{p}$
and $\lambda_{p}$ to the Cheeger problem in the limit $p\rightarrow 1$.

######  Theorem 2.4.

Let $M\subseteq\mathbb{R}^{d}$ be a bounded manifold of dimension $d$ with
Lipschitz boundary. Then in the limit $p\rightarrow 1$

  1. (a)

the first eigenvalue $\lambda_{p}$ converges to $h(M)$;

  2. (b)

after possibly passing to a subsequence, the first eigenfunction of
$-\Delta_{p}$ converges in $L^{1}$ to a function $u\in BV(M)$ (the space of
functions of bounded variation on $M$) such that the superlevel sets
$A_{t}:=\\{x\in M\mid u(x)>t\\}$ are either null sets or Cheeger sets of $M$
for almost all $t>0$;

  3. (c)

if $M$ has a unique Cheeger set $A$, then $u_{p}$ converges in $L^{1}$ to a
scalar multiple of the characteristic function $\chi_{A}$.

###### Proof sketch..

For (a) see e.g. [20, Corollary 6]. For (b) one can follow [20, Theorem 8,
Remark 10] to show that by compactness of the embedding $BV(M)\rightarrow
L^{1}(M)$ there must be a subsequence converging to some $u$ in $L^{1}$, for
which one can deduce from (a) that it must fulfil

$\int_{0}^{\infty}\ell_{d-1}(\partial
A_{t})\;dt=\int_{0}^{\infty}h(M)\ell_{d}(A_{t})\;dt.$ (14)

As $\ell_{d-1}(\partial A_{t})\leq h(M)\ell_{d}(A_{t})$, this implies that
$\ell_{d-1}(\partial A_{t})=h(M)\ell_{d}(A_{t})$ for almost all $t$, and hence
$A_{t}$ is either a null set or a Cheeger set. For (c) we refer to [29, Remark
3.2.], where it is noted that one can show $\ell_{d-1}(\partial
A_{t})=h(M)\ell_{d}(A)$ for _all_ $t$. By uniqueness of the Cheeger set, this
means that $A_{t}=A$ for all $t>0$ where $A_{t}$ is not a null set and that
this implies that $u$ is a suitably scaled characteristic function on a
Cheeger set of $M$. ∎

In section 3.2 we will show an analogous statement to 2.4(a) for the dynamic
$p$-Laplacian introduced in section 3. 2.4(b) and (c) are investigated
numerically in section 5, in particular we study the improvement of the
Cheeger ratios of the level sets as $p\to 1$.

##  3 A dynamic $p$-Laplacian

###  3.1 Motivation

We now aim to generalise the $p$-Laplacian to a dynamic $p$-Laplacian which
shares the same connections to the dynamic Cheeger problem as the classic
$p$-Laplacian does to the classic Cheeger problem. Recall that we consider a
volume-preserving diffeomorphism $T:M\rightarrow M$, use the pushforward
operator $T_{*}:L^{p}(M)\rightarrow L^{p}(M)$, $T_{*}u:=u\circ T^{-1}$, and
its adjoint $T^{*}:L^{q}(M)\rightarrow L^{q}(M)$, $T^{*}u:=u\circ T$, the
Koopman operator. Note that for $u\in W^{1,p}_{0}(M)$ we have $T_{*}u\in
W^{1,p}_{0}(M)$ by compactness of $M$ so we can use $T_{*}$ as an operator
from $W^{1,p}_{0}(M)$ to $W^{1,p}_{0}(M)$.

Equipped with this we define a dynamic version
${F}^{D}:W^{1,p}_{0}(M)\rightarrow\mathbb{R}$ of $F$ as

$\displaystyle{F}^{D}(u)$ $\displaystyle:=\frac{1}{2}\left(\|\nabla
u\|_{p}^{p}+\|\nabla T_{*}u\|_{p}^{p}\right)$ (15)

or equivalently

${F}^{D}=\frac{1}{2}(F+F\circ T_{*}).$

It now follows directly from the differentiability of $F$ and boundedness of
$T_{*}$ that ${F}^{D}$ is differentiable with

$\displaystyle({F}^{D})^{\prime}(u)v$
$\displaystyle=\frac{1}{2}\bigl{(}F^{\prime}(u)v+F^{\prime}(T_{*}u)T_{*}v\bigr{)}$
(16)
$\displaystyle=\frac{1}{2}p\bigl{(}(\Delta_{p}u)v+(T^{*}\Delta_{p}T_{*}u)v\bigr{)}.$
(17)

To get from (16) to (17) we see the Koopman operator $T^{*}$ as the purely
linear algebraic adjoint of $T_{*}$, meaning that $T^{*}$ maps some
$f\in(W^{1,p}_{0}(M))^{*}$ to $f\circ T_{*}$. This directly implies
$F^{\prime}(T_{*}u)T_{*}v=(F^{\prime}(T_{*}u)\circ
T_{*})v=(T^{*}F^{\prime}(T_{*}u))v$. While this differs from the common
definition of $T^{*}$ as a map on $L^{q}(M)$ defined by $f\mapsto f\circ T$,
the two definitions coincide on $L^{q}$ under the usual identification of some
$f\in L^{q}$ with the functional
$(u\mapsto\int_{M}uf)\in(W^{1,p}_{0}(M))^{*}$, as by volume preservation of
$T$:

$\left(u\mapsto\int_{M}u\cdot f\right)\circ
T_{*}=\left(u\mapsto\int_{M}(u\circ T^{-1})\cdot
f\right)=\left(u\mapsto\int_{M}u\cdot(f\circ T)\right).$

This motivates the following definition for a dynamic $p$-Laplacian:

$\displaystyle{\Delta}^{D}_{p}u$
$\displaystyle:=-\frac{1}{p}({F}^{D})^{\prime}(u)=\frac{1}{2}\left(\Delta_{p}u+T^{*}\Delta_{p}T_{*}u\right).$
(18)

Plugging in $F^{\prime}(u)v=p\int_{M}|\nabla u|^{p-2}\nabla u\nabla v$ into
(16) we get the equivalent definition

$({\Delta}^{D}_{p}u)v:=\frac{1}{2}\left(\int_{M}\nabla u|\nabla u|^{p-2}\nabla
v+\int_{M}\nabla T_{*}u|\nabla T_{*}u|^{p-2}\nabla T_{*}v\right).$ (19)

Analogously to 2.1, we have:

######  Theorem 3.1.

Let $d\geq 2$, $1<p<2$ and $M\subseteq\mathbb{R}^{d}$ a compact
$d$-dimensional manifold with nonempty Lipschitz boundary. Define
${J}^{D}:W_{0}^{1,p}(M)\backslash\\{0\\}\rightarrow\mathbb{R}$ by

${J}^{D}(u):=\frac{{F}^{D}(u)}{G(u)}.$

The infimum

${\lambda}^{D}_{p}(M,T):=\inf_{u\in
W_{0}^{1,p}(M)\backslash\\{0\\}}{J}^{D}(u)$ (20)

is positive and achieved by a function ${u}^{D}_{p}\in W_{0}^{1,p}(M)$
satisfying

$-{\Delta}^{D}_{p}{u}^{D}_{p}={\lambda}^{D}_{p}{u}^{D}_{p}|{u}^{D}_{p}|^{p-2}$
(21)

weakly. We call ${\lambda}^{D}_{p}$ the _first eigenvalue_ of
$-{\Delta}^{D}_{p}$ and ${u}^{D}_{p}$ a _first eigenfunction_ of
$-{\Delta}^{D}_{p}$.

######  Remark 3.1.

We do not prove a uniqueness result as in 2.3 and thus do not talk about _the_
first eigenfunction of $-{\Delta}^{D}_{p}$ in 3.1. The eigenfunction
$u^{D}_{p}$ is also not proven to be positive, in analogy to the positivity
result in 2.2.

###### Proof of 3.1.

We use the standard direct method that is also used for the (static)
$p$-Laplacian. The crucial part is weak lower semicontinuity of the numerator
of ${J}^{D}$. It is well known that $u\mapsto\|\nabla u\|_{p}^{p}$ is weakly
lower semicontinuous, as it is continuous and convex (see e.g. [6, Theorem
1.5.3]). To establish weak lower semicontinuity of $u\mapsto\|\nabla
T_{*}u\|_{p}^{p}$, we first show that for every weakly convergent series
$u_{k}\rightharpoonup u$ we also have $T_{*}u_{k}\rightharpoonup T_{*}u$. This
is a simple consequence of linearity and boundedness of $T_{*}$: For any $f\in
W^{1,p}_{0}(M)^{*}$ we can define $\tilde{f}:=f\circ T_{*}\in
W^{1,p}_{0}(M)^{*}$, and hence $u_{k}\rightharpoonup u$ implies

$f(T_{*}u_{k})=\tilde{f}(u_{k})\rightarrow\tilde{f}(u)=f(T_{*}u).$

This proves $T_{*}u_{k}\rightharpoonup T_{*}u$ and implies that for a weakly
lower semicontinuous map $H$, the map $H\circ T_{*}$ is also weakly lower
semicontinuous. In particular the functional $u\mapsto\|\nabla
T_{*}u\|_{p}^{p}$ and hence also the functional

${F}^{D}(u):=\frac{1}{2}\left(\|\nabla u\|_{p}^{p}+\|\nabla
T_{*}u\|_{p}^{p}\right)$

are weakly lower semicontinuous. Now we can proceed with a standard argument,
as seen e.g. in [6, Theorem 2.6.11]: Clearly, ${J}^{D}$ is bounded from below
by $0$ and thus has a non-negative finite infimum ${\lambda}^{D}_{p}$. The
Poincaré inequality establishes that for some positive $C$,

$G(u)=\|u\|_{p}^{p}\leq C\|\nabla u\|_{p}^{p}\leq 2C{F}^{D}(u)$

and thus ${\lambda}^{D}_{p}\geq\frac{1}{2C}$ is positive. Let $(u_{k})_{k}$ be
a minimising sequence for ${J}^{D}$. We can assume without loss that
$\|u_{k}\|_{p}=1$ for all $k$, as ${J}^{D}$ is homogeneous. The sequence is
therefore bounded in $W_{0}^{1,p}(M)$ and we can pass to a weakly convergent
subsequence and assume $u_{k}\rightharpoonup u$ for some $u\in
W_{0}^{1,p}(M)$. Now by the Rellich–Kondrachov theorem (see e.g. [2, Theorem
6.3]) the inclusion mapping $W_{0}^{1,p}(M)\xhookrightarrow{}L^{p}(M)$ is
compact, and as compact operators map weakly convergent sequences to strongly
convergent ones, we have $u_{k}\rightarrow u$ in $L^{p}$ and hence
$G(u):=\|u\|_{p}^{p}=1$. Finally, by weak lower semicontinuity of ${F}^{D}$ we
have

${\lambda}^{D}_{p}\leq{J}^{D}(u)={F}^{D}(u)\leq\liminf_{k}{F}^{D}(u_{k})=\liminf_{k}{J}^{D}(u_{k})={\lambda}^{D}_{p}$

and hence ${\lambda}^{D}_{p}$ is achieved by some ${u}^{D}_{p}:=u$. To show
that ${u}^{D}_{p}$ satisfies (21), note that since ${u}^{D}_{p}$ is a global
minimum of ${J}^{D}={F}^{D}/G$ we have $({J}^{D})^{\prime}({u}^{D}_{p})=0$, so
for all $v\in W_{0}^{1,p}(M)$:

$\displaystyle 0$ $\displaystyle=({J}^{D})^{\prime}({u}^{D}_{p})(v)$
$\displaystyle=\frac{1}{G({u}^{D}_{p})^{2}}(G({u}^{D}_{p})({F}^{D})^{\prime}({u}^{D}_{p})(v)-{F}^{D}({u}^{D}_{p})G^{\prime}({u}^{D}_{p})(v))$
$\displaystyle=({F}^{D})^{\prime}({u}^{D}_{p})(v)-{\lambda}^{D}_{p}G^{\prime}({u}^{D}_{p})(v).$

Now the last expression is the weak form of
$-p{\Delta}^{D}_{p}{u}^{D}_{p}=p{\lambda}^{D}_{p}{u}^{D}_{p}|{u}^{D}_{p}|^{p-2}$,
which shows the claim. ∎

###  3.2 Behaviour as $p\rightarrow 1$

2.4(a) states that the first eigenvalue $\lambda_{p}$ of $-\Delta_{p}$ with
Dirichlet boundary conditions converges to $h(M)$ as $p\rightarrow 1$ in the
static situation. We prove the analogous statement for the _dynamic_
$p$-Laplacian with Dirichlet boundary conditions. In the proof, an inequality
that plays a crucial role is

$\lambda_{p}\geq\left(\frac{h(M)}{p}\right)^{p},$

which is called the _Cheeger_ inequality in the case $p=2$ (see [10] for $p=2$
and [21, Appendix] or [20, Theorem 3] for general $p$). We prove this
inequality for the _dynamic_ $p$-Laplacian with Dirichlet boundary conditions.
For the readers’ convenience, we state two classical results that we use in
the proof.

######  Theorem 3.2 (Co-area formula).

Let $M\subseteq{\mathbb{R^{d}}}$ be a compact $d$-dimensional submanifold with
Lipschitz boundary and $u\in C_{0}^{\infty}(M)$. Then

$\int_{M}|\nabla u|=\int_{-\infty}^{\infty}\ell_{d-1}(\partial A(t))\ dt,$

where $A(t)=\\{x\in M\mid u(x)>t\\}$ and the boundary is relative to $M$.

###### Proof.

See e.g. [9, Theorem I.3.4 and Corollary I.3.1]. The formula is commonly
stated with $u^{-1}(t)$ instead of $\partial A(t)$. For a proof that
$u^{-1}(t)=\partial A(t)$ for almost all $t$ see A.1. ∎

The second theorem is a special case of Cavalieri’s Principle:

######  Theorem 3.3 (Cavalieri’s Principle).

Let $M\subseteq{\mathbb{R}^{d}}$ be a compact $d$-dimensional manifold with
Lipschitz boundary and $f:M\rightarrow\mathbb{R}$ be a non-negative measurable
function. Then

$\int_{M}f=\int_{0}^{\infty}\ell_{d}(\\{f>t\\})\ dt.$

###### Proof.

See [9, Proposition I.3.3] ∎

Equipped with these results, we now proceed to the proof:

######  Theorem 3.4.

Let $M\subset\mathbb{R}^{d}$ be a compact $d$-dimensional submanifold with
Lipschitz boundary, $T:M\rightarrow M$ a smooth volume-preserving
diffeomorphism and $p\geq 1$. Then

${\lambda}^{D}_{p}(M,T)\geq\left(\frac{{h}^{D}(M,T)}{p}\right)^{p}.$ (22)

###### Proof.

We modify the proof that is given in the appendix of [21] and [20, Theorem 3].

Step 1. First we prove the inequality for $p=1$. For this, assume $w\in
C^{\infty}_{0}(M)$ and define $A(t):=\\{x\in M\mid w(x)>t\\}$. By the coarea
formula

$\int_{M}|\nabla w|=\int_{-\infty}^{\infty}\ell_{d-1}(\partial A(t))\ dt,$
(23)

where the boundary is relative to $M$. If we apply the coarea formula to
$T_{*}w$ we get

$\displaystyle\int_{M}|\nabla T_{*}w|$
$\displaystyle=\int_{-\infty}^{\infty}\ell_{d-1}(\partial\\{x\in M\mid
T_{*}w(x)>t\\})\ dt.$ (24)
$\displaystyle=\int_{-\infty}^{\infty}\ell_{d-1}(\partial\\{x\in M\mid
w(T^{-1}(x))>t\\})\ dt.$ (25)
$\displaystyle=\int_{-\infty}^{\infty}\ell_{d-1}(\partial\\{Ty\mid w(y)>t\\})\
dt.$ (26) $\displaystyle=\int_{-\infty}^{\infty}\ell_{d-1}(\partial\\{Ty\mid
y\in A(t)\\})\ dt.$ (27)
$\displaystyle=\int_{-\infty}^{\infty}\ell_{d-1}(\partial(TA(t)))\ dt.$ (28)

Now by A.1(c) we can express the integrals on the right-hand side of (23) and
(28) by using the sets $B(t):=\\{x\in M\mid|w(x)|>t\\}$ to get

$\displaystyle\int_{M}|\nabla w|$
$\displaystyle=\int_{0}^{\infty}\ell_{d-1}(\partial B(t))\ dt.$ (29)
$\displaystyle\int_{M}|\nabla T_{*}w|$
$\displaystyle=\int_{0}^{\infty}\ell_{d-1}(\partial(TB(t)))\ dt.$ (30)

Combining (29) and (30) we can continue with

$\displaystyle\int_{M}\frac{1}{2}(|\nabla w|+|\nabla T_{*}w|)$
$\displaystyle=\frac{1}{2}\int_{0}^{\infty}\ell_{d-1}(\partial
B(t))+\ell_{d-1}(\partial TB(t))\ dt$ (31)
$\displaystyle=\frac{1}{2}\int_{0}^{\max|w|}\frac{\ell_{d-1}(\partial
B(t))+\ell_{d-1}(\partial TB(t))}{\ell_{d}(B(t))}\ell_{d}(B(t))dt$ (32)
$\displaystyle\geq{h}^{D}(M,T)\int_{0}^{\max|w|}\ell_{d}(B(t))\ dt$ (33)
$\displaystyle={h}^{D}(M,T)\int_{M}|w|,$ (34)

where in the last step we have applied Cavalieri’s principle to $|w|$.

For the step (32) to (33) we also need to justify that $B(t)$ is in the
admissible class for the Cheeger problem for almost all $t$. First observe
that, as $|w|$ is continuous, and vanishes on the boundary of $M$, $B(t)$ does
not touch $\partial M$ if $t>0$. As for the regularity of $\partial B(t)$,
note that by A.1(b) we have $\partial B(t)$ is the disjoint union of $\partial
A(t)$ and $\partial A(-t)$ for almost all $t$. Sard’s theorem tells us that
$w^{-1}(t)$ contains no critical points for almost all $t$. Hence, by the
implicit function theorem, $A(t)$ has $C^{\infty}$ boundary for almost all $t$
and so does $B(t)$. After applying the same argument to $TB(t)$ we know that
we may indeed bound the fraction in (32) by ${h}^{D}(M,T)$.

The inequality (34) also holds for $w\in W_{0}^{1,1}(M)$ as
$C^{\infty}_{0}(M)$ is dense in $W^{1,1}_{0}(M)$. This shows the claim for
$p=1$.

Step 2. Next we generalise to $p>1$, essentially by reducing the result to the
case in Step 1. For $p>1$ and $v\in W^{1,p}(M)$ define $\Phi(v):=|v|^{p-1}v$.
As $(|x|^{p-1}x)^{\prime}=p|x|^{p-1}$, we can use the chain rule to get that
the gradient of $\Phi(v)$ is

$\nabla\Phi(v(x))=p|v(x)|^{p-1}\nabla v(x).$

By Hölder’s inequality we can estimate

$\displaystyle\int_{M}|\nabla\Phi(v)|=p\int_{M}|v|^{p-1}|\nabla v|\leq
p\left(\int_{M}|v|^{(p-1)\frac{p}{p-1}}\right)^{\frac{p-1}{p}}\|\nabla
v\|_{p}=p\|v\|_{p}^{p-1}\|\nabla v\|_{p},$ (35)

implying that $w:=\Phi(v)\in W^{1,1}_{0}(M)$ for $v\in W^{1,p}_{0}(M)$.

We repeat for $T_{*}v$. First note that

$\Phi(T_{*}v(x))=\Phi(v(T^{-1}(x)))=|v(T^{-1}x)|^{p-1}v(T^{-1}x)=\left[T_{*}(|v|^{p-1}v)\right](x)=T_{*}(\Phi(v))(x)$

and hence $\Phi(T_{*}v)=T_{*}(\Phi(v))$. Thus, as above:

$\int_{M}|\nabla T_{*}w|=\int_{M}|\nabla
T_{*}\Phi(v)|=\int_{M}|\nabla\Phi(T_{*}(v))|\leq p\|T_{*}v\|^{p-1}_{p}\|\nabla
T_{*}v\|_{p}=p\|v\|^{p-1}_{p}\|\nabla T_{*}v\|_{p}$ (36)

using that $T$ is volume-preserving in the last step. Now on the right-hand
side of (36), both norms are bounded by compactness of $M$ and smoothness of
$T$. Hence, also $T_{*}w\in W^{1,1}_{0}(M)$. We can therefore apply inequality
(34) to get

$\displaystyle{h}^{D}(M,T)$
$\displaystyle\overset{\eqref{coareacavalieri}}{\leq}\frac{\int_{M}|\nabla
w|+|\nabla T_{*}w|}{2\int_{M}|w|}$ (37)
$\displaystyle\overset{\eqref{holder},\eqref{holder2}}{\leq}\frac{p\|v\|^{p-1}_{p}(\|\nabla
v\|_{p}+\|\nabla T_{*}v\|_{p})}{2\int_{M}|v|^{p}}$ (38)
$\displaystyle=p\frac{\|\nabla v\|_{p}+\|\nabla T_{*}v\|_{p}}{2\|v\|_{p}}$
(39)

Raising this to the power of $p$ and using that convexity of $x^{p}$ implies
$(\frac{a+b}{2})^{p}\leq\frac{a^{p}+b^{p}}{2}$ we arrive at

$\displaystyle\left(\frac{{h}^{D}(M,T)}{p}\right)^{p}$
$\displaystyle\leq\left(\frac{\|\nabla v\|_{p}+\|\nabla
T_{*}v\|}{2\|v\|_{p}}\right)^{p}$ (40) $\displaystyle\leq\frac{\|\nabla
v\|_{p}^{p}+\|\nabla T_{*}v\|_{p}^{p}}{2\|v\|^{p}_{p}}$ (41)

which gives the desired inequality after passing to the infimum. ∎

Finally, we can turn our attention to the behaviour of ${\lambda}^{D}_{p}$ as
$p\to 1$.

######  Theorem 3.5.

Let $M\subseteq\mathbb{R}^{d}$ and $T:M\rightarrow M$ be as in 3.4. Then the
first eigenvalue ${\lambda}^{D}_{p}(M,T)$ converges to ${h}^{D}(M,T)$ as
$p\rightarrow 1$.

###### Proof.

We use a mixture of arguments from [20, Corollary 6], and [14, Lemma A.1 and
Theorem 3.1]. By Equation 22 we already have a lower bound on
${\lambda}^{D}_{p}(M,T)$ that converges to ${h}^{D}(M,T)$ as $p\rightarrow 1$.
To construct an upper bound, consider a subdomain $A_{k}\subset M$ with smooth
boundary not touching $\partial M$ and

$\frac{\ell_{d-1}(\partial A_{k})+\ell_{d-1}(\partial
TA_{k})}{2\ell_{d}(A_{k})}-{h}^{D}(M,T)<\frac{1}{k}\quad\text{for a
$k\in\mathbb{N}$}.$

Define

$\displaystyle f_{\varepsilon}(x):=\begin{cases}1&x\in A_{k}\\\
1-\frac{1}{\varepsilon}d(x,A_{k})&0<d(x,A_{k})<\varepsilon\\\
0&d(x,A_{k})\geq\varepsilon\\\ \end{cases}$ (42)

with $d(x,A_{k}):=\inf_{y\in A_{k}}|x-y|$. Then $f_{\varepsilon}$ is in
$W_{0}^{1,p}(M)$ for $\varepsilon$ small enough and

${\lambda_{p}}^{D}(M,T)\leq\frac{\|\nabla f_{\varepsilon}\|_{p}^{p}+\|\nabla
T_{*}f_{\varepsilon}\|^{p}_{p}}{2\|f_{\varepsilon}\|_{p}^{p}}=:C_{p,\varepsilon,k}.$
(43)

To remove the dependence on $p$, first note that $|\nabla f_{\varepsilon}|\leq
G_{\varepsilon}$ for some constant $G_{\varepsilon}$ and $|\nabla
f_{\varepsilon}|^{p}$ converges pointwise to $|\nabla f_{\varepsilon}|$ for
$p\rightarrow 1$, hence by compactness of $M$ and dominated convergence
$\lim_{p\to 1}\|\nabla f_{\varepsilon}\|_{p}^{p}=\|\nabla
f_{\varepsilon}\|_{1}$. Using that $T$ is continuously differentiable and
applying the same argument to $\|\nabla T_{*}f_{\varepsilon}\|_{p}^{p}$ and
$\|f_{\varepsilon}\|_{p}^{p}$ yields

$\lim_{p\rightarrow 1}C_{p,\varepsilon,k}=\frac{\|\nabla
f_{\varepsilon}\|_{1}+\|\nabla
T_{*}f_{\varepsilon}\|_{1}}{2\|f_{\varepsilon}\|_{1}}=:C_{\varepsilon,k}.$

Now $\lim_{\varepsilon\rightarrow 0}\|f_{\varepsilon}\|_{1}=\ell_{d}(A_{k})$
by dominated convergence and $\lim_{\varepsilon\rightarrow 0}\|\nabla
f_{\varepsilon}\|_{1}=\ell_{d-1}(\partial A_{k})$, as $|\nabla
f_{\varepsilon}|=\frac{1}{\varepsilon}$ almost everywhere on
$\Gamma_{\varepsilon}:=\\{x\in M\backslash A\mid d(x,M)\leq\varepsilon\\}$ and
zero elsewhere so $\lim_{\varepsilon\rightarrow 0}\int_{M}|\nabla
f_{\varepsilon}|=\lim_{\varepsilon\rightarrow
0}\ell_{d}(\Gamma_{\varepsilon})/\varepsilon=\ell_{d-1}(\partial A)$ (see also
proof of [14, Theorem 3.1]). Following the proof of [14, Theorem 3.1], using
[14, Lemma A.1] we have $\lim_{\varepsilon\rightarrow 0}\|\nabla
T_{*}f_{\varepsilon}\|_{1}=\ell_{d-1}(T\partial A_{k})=\ell_{d-1}(\partial
TA_{k})$, which leads us to

$\lim_{\varepsilon\rightarrow 0}C_{\varepsilon,k}=\frac{\ell_{d-1}(\partial
A_{k})+\ell_{d-1}(\partial TA_{k})}{2\ell_{d}(A_{k})}=:C_{k}.$

Finally by definition of $A_{k}$ we have

$\lim_{k\rightarrow\infty}C_{k}={h}^{D}(M,T).$ (44)

Now we can conclude by taking successive limits. By (43) we have for fixed $k$
and $\varepsilon$

$\limsup_{p\rightarrow 1}{\lambda}^{D}_{p}(M,T)\leq C_{\varepsilon,k},$

which, after passing to the limit $\varepsilon\rightarrow 0$ implies that for
all $k$

$\limsup_{p\rightarrow 1}{\lambda}^{D}_{p}(M,T)\leq C_{k}.$

A final limit $k\rightarrow\infty$ together with Equation 44 shows

$\limsup_{p\rightarrow 1}{\lambda}^{D}_{p}(M,T)\leq{h}^{D}(M,T).$

Recalling that Equation 22 gave us $\liminf_{p\rightarrow
1}\lambda_{p}(M,T)\geq{h}^{D}(M,T)$ we can conclude that

$\lim_{p\rightarrow 1}{\lambda}^{D}_{p}(M,T)={h}^{D}(M,T).$

∎

##  4 Numerical approximation of eigenfunctions

Finding eigenfunctions of the $p$-Laplacian computationally is much harder for
$p\neq 2$ than for $p=2$. In the latter case, the operator is the standard
linear Laplacian and finite element methods can readily be applied [34]. For
the $p$-Laplacian, Horák [19] uses a modified gradient descent for the first
eigenfunction and a mountain pass algorithm for the second one. Lefton and Wei
[21] use a penalty formulation together with the Levenberg–Marquard algorithm
to approximate the first eigenfunction. For our purpose, we will modify a
method by Yao and Zhou [36], as it is general enough to be easily applied to
our dynamical functional ${F}^{D}$.

###  4.1 Abstract algorithm for the first eigenpair

We will use the algorithm mostly for the first eigenfunction. As it simplifies
substantially in this case and is essentially a descent algorithm, we
summarise this special case here. Recall that for the first eigenpair
$(u_{p},\lambda_{p})$ of the $p$-Laplacian, we have

$\lambda_{p}=\inf_{u\in W^{1,p}_{0}(M)\backslash\\{0\\}}J(u)$

and the infimum is attained at the first eigenfunction $u_{p}$. For the
descent algorithm, [36] proposes to use the descent direction $-d\in
W^{1,q}_{0}(M)$ that fulfils

$\int_{M}\nabla d\cdot\nabla v=J^{\prime}(u)(v)$ (45)

for all $v\in W^{1,p}_{0}(M)$. We will denote this descent direction by
$\operatorname{grad}J(u):=d$. In the case of $C^{1}$ boundary, it is shown in
[33, Theorem 7.5] that such a $d$ exists.444To apply [33, Thm. 7.5], set
$B[\mathbf{u},\Phi]:=\int_{M}\nabla\mathbf{u}\nabla\Phi$,
$F(\Phi):=J^{\prime}(u)(\Phi)$ (using the respective notations $B,\mathbf{u}$
and $\Phi$ from [33] and $u,F,J^{\prime}$ from our definitions) and note that
the set $N_{q}$ from the assumptions of [33, Thm. 7.5.] equals $\\{0\\}$,
hence trivially $F(z)=0$ for all $z\in N_{q}$. Note that because $M$ is
bounded and $p<q$ we have $W^{1,q}_{0}(M)\subseteq W^{1,2}_{0}(M)\subseteq
W^{1,p}_{0}(M)$, so we also get $d\in W^{1,p}_{0}(M)$. Note further that
$J^{\prime}(u)(d)=\|\nabla d\|_{2}^{2}$ and thus

$\displaystyle J(u-td)$ $\displaystyle=J(u)-t\,J^{\prime}(u)(d)+o(t)$ (46)
$\displaystyle=J(u)-t\|\nabla d\|_{2}^{2}+o(t)$ (47)
$\displaystyle<J(u)-t\|\nabla d\|_{2}^{2}/4$ (48)

for small $t$, i.e. $-d$ is a descent direction and a suitable stepsize $t$
can be found by an Armijo-type stepsize search.

The corresponding descent scheme is given in Algorithm 1, which is based on
the local min-max algorithm from [36, Section 3], specialised for the first
eigenfunction and with the modifications from [36, Section 4.1] filled in.

Algorithm 1 Computation of the first eigenfunction

  1. 1.

Set $k=1$, choose an initial guess $u^{[1]}\in W^{1,p}_{0}(M)$ and a parameter
$\alpha>0$.

  2. 2.

While $\|\operatorname{grad}J(u^{[k]})\|_{p}$ is too large do

    1. (a)

Compute the descent direction $d^{[k]}:=-\operatorname{grad}J(u^{[k]})$.

    2. (b)

For $s\in\mathbb{R}$ define

$u^{[k+1]}(s):=\frac{u^{[k]}+sd^{[k]}}{\|u^{[k]}+sd^{[k]}\|_{p}}.$

    3. (c)

Set $s^{*}:=\alpha/\max(1,\|d^{[k]}\|_{p})$ and check whether

$J(u^{[k+1]}(s^{*}))\leq J(u^{[k]})-\frac{s^{*}}{4}\|\nabla d^{[k]}\|_{2}^{2}$

is fulfilled. If it is not, keep halving $s^{*}$ until it is.

    4. (d)

Set $u^{[k+1]}:=u^{[k+1]}(s^{*})$.

    5. (e)

Increment $k\to k+1$.

    6. (f)

If $s^{*}$ is smaller than some tolerance, go to Step 3.

  3. 3.

Return $u_{p}:=u^{[k]}$.

###  4.2 Application to the static and the dynamic $p$-Laplacian

#### The functionals.

For the static $p$-Laplacian $\Delta_{p}$ we have $F(u)=\|\nabla u\|_{p}^{p}$
and $G(u)=\|u\|_{p}^{p}$ so that the derivative of $J$ is

$\displaystyle J^{\prime}(u)(v)$
$\displaystyle=\frac{1}{G(u)^{2}}(G(u)F^{\prime}(u)(v)-F(u)G^{\prime}(u)(v))$
(49) $\displaystyle=\frac{p}{b^{2}}\int_{M}(b|\nabla u|^{p-2}\nabla u\nabla
v-a|u|^{p-2}uv),$ (50)

where $a=\|\nabla u\|_{p}^{p}$ and $b=\|u\|_{p}^{p}$ [36, Section 4.1]. Hence,
the descent direction $d$ is given by the solution of

$\int_{M}\nabla d\cdot\nabla v=\frac{p}{b^{2}}\int_{M}(b|\nabla u|^{p-2}\nabla
u\nabla v-au|u|^{p-2}v)\qquad\text{for all }v\in W_{0}^{1,p}(M).$

In case of the dynamic $p$-Laplacian ${\Delta}^{D}_{p}$, we consider

${J}^{D}(u):=\frac{{F}^{D}(u)}{G(u)}=\frac{1}{2}\frac{F(u)+F(T_{*}u)}{G(u)}.$

As in (16) we have

$({F}^{D})^{\prime}(u)v=\frac{1}{2}(F^{\prime}(u)v+F^{\prime}(T_{*}u)T_{*}v),$

so the derivative of ${J}^{D}$ is

$\displaystyle({J}^{D})^{\prime}(u)v$
$\displaystyle=\frac{1}{G(u)^{2}}(G(u)({F}^{D})^{\prime}(u)v-{F}^{D}(u)G^{\prime}(u)v)$
(51)
$\displaystyle=\frac{1}{b^{2}}\left(\frac{1}{2}b(F(u)v+F(T_{*}u)T_{*}v)-{a}^{D}G^{\prime}(u)v\right)$
(52) $\displaystyle=\frac{p}{b^{2}}\int_{M}\left(\frac{1}{2}\left(|\nabla
u|^{p-2}\nabla u\cdot\nabla v+|\nabla(u\circ T)|^{p-2}\nabla(u\circ
T)\cdot\nabla(v\circ T)\right)-{a}^{D}|u|^{p-2}uv\right),$ (53)

where ${a}^{D}:=\frac{1}{p}{F}^{D}(u)$. Defining $A(x):=DT^{-T}(x)$ and using
volume-preservation of $T$ we make the replacement

$\int_{M}|\nabla(u\circ T)|^{p-2}\nabla(u\circ T)\cdot\nabla(v\circ
T)=\int_{M}|A\nabla u|^{p-2}(A\nabla u)\cdot(A\nabla v).$

This yields

$({J}^{D})^{\prime}(u)v=\frac{p}{b^{2}}\int_{M}\left(\frac{1}{2}\left(|\nabla
u|^{p-2}\nabla u\cdot\nabla v+|A\nabla u|^{p-2}(A\nabla u)\cdot(A\nabla
v)\right)-{a}^{D}|u|^{p-2}uv\right).$

#### The descent direction.

We found that in the case of the dynamic $p$-Laplacian, it is not advisable to
define $d$ using (45) directly replacing $J^{\prime}$ with
$({J}^{D})^{\prime}$, as the resulting $d$ exhibits high irregularity in areas
where the norm of the derivative of $T$ is large, leading to undesirably small
steps in the descent algorithm, similar to the problems gradient descent
exhibits when not conditioned properly.

As a remedy, we instead use ${d}^{D}$ defined by

$\int_{M}\nabla{d}^{D}\cdot\left(\frac{1}{2}(I+A^{T}A)\nabla
v\right)=({J}^{D})^{\prime}(u)v\qquad\text{for all }v\in W^{1,p}_{0}(M)$ (54)

and denote this choice by ${\operatorname{grad}}^{D}{J}^{D}:={d}^{D}$. Note
that because

$\Big{\langle}\ell,\frac{1}{2}(I+A(x)^{T}A(x))\ell\Big{\rangle}\geq\frac{1}{2}|\ell|^{2}$

for all $x\in M$ and $\ell\in\mathbb{R}^{d}$, and the coefficients of $A$ are
bounded in $M$, the left-hand side of (54) defines a uniformly strongly
elliptic regular Dirichlet bilinear form of order $1$, so [33, Theorem 7.5] is
again applicable if $M$ has $C^{1}$ boundary.

In our experiments, this choice of ${d}^{D}$ lead to higher regularity of the
descent direction, resulting in considerably bigger step sizes (see Figure 3
for a visual comparison of $\operatorname{grad}{J}^{D}$ and
${\operatorname{grad}}^{D}{J}^{D}$).

Figure 3: The gradients $\operatorname{grad}{J}^{D}(u)$ (left) and
${\operatorname{grad}}^{D}{J}^{D}(u)$ (right) at the starting point $u$ of
Algorithm 1 for $p=1.83$ in the example in section 5.2. As a starting point
$u$ we choose the first eigenfunction of ${\Delta}^{D}_{2}$. Note the log
transformation in the left plot, which was used in order to visualise the
different orders of magnitude.

#### Initialization.

Finally, one needs to find an initial guess for $u^{[1]}$ in the first step of
Algorithm 1. Yao and Zhou suggest finding functions with appropriate nodal
line structures (see [36, Remark 3.1(a)]). We used the first eigenfunction of
${\Delta}^{D}_{2}$ for this purpose, which we computed using the Cauchy–Green
approach from [16, Section 3.1].

###  4.3 The second eigenpair

One of our numerical examples in section 5 is on the torus $\mathbb{T}^{2}$,
which has no boundary. This has the consequence that the infimum of ${J}^{D}$
is $0$, realised by any constant function and carrying no information about
the dynamics. We thus need to compute the second eigenfunction and therefore
briefly describe the algorithm of Yao and Zhou for this case and present our
version for the dynamic $p$-Laplacian $\smash{{\Delta}^{D}_{p}}$.

Note that the first eigenvalue $\lambda_{p}=\lambda_{p}^{(1)}$ of $\Delta_{p}$
is isolated [24, Theorem 9]. As the spectrum of $\Delta_{p}$ is closed [24,
Theorem 3], we may speak of the second eigenvalue $\lambda^{(2)}_{p}$ of
$\Delta_{p}$ as the infimum over all eigenvalues bigger than
$\lambda_{p}^{(1)}$. We assume this to be true for the dynamic $p$-Laplacian
$\smash{{\Delta}^{D}_{p}}$ as well. The second eigenvalue can be characterised
by the min-max-principle [36, Equation (2.1)](cf. Appendix A.1 for a
derivation in the linear case)

$\lambda^{(2)}_{p}=\min_{v\in L^{\prime}}\max_{u\in[L,v]\cap S}J(u),$ (55)

where $L$ is the span of the first eigenfunction $u_{p}=u_{p}^{(1)}$,
$L^{\prime}$ is a complement of $L$ (i.e. such that $W^{1,p}_{0}(M)=L\oplus
L^{\prime}$) and $[L,v]:=\textrm{span}(L,v)$. By $S:=\\{u\in
W^{1,p}_{0}(M)\mid\|u\|=1\\}$, we denote the unit sphere in $W^{1,p}_{0}(M)$
and define $[L,v]_{S}:=\\{t_{1}u_{p}^{(1)}+t_{2}v\mid
t_{1}^{2}+t_{2}^{2}=1\\}$. Since $J(tu)=u$ for all $t\neq 0$, and
$[L,v]_{S}\xrightarrow{{u\ \mapsto\ u/\|u\|}}[L,v]\cap S$ is a bijection, we
have

$\max\\{J(u)\mid u\in[L,v]_{S}\\}=\max\\{J(u)\mid u\in[L,v]\cap S\\}.$ (56)

We work directly with the left-hand side of (56) as the optimisation is now
finite-dimensional over $t_{1},t_{2}$. From a maximiser of the left-hand side
of (56) we can immediately obtain a maximiser of the right-hand side.

A local optimum of the maximisation in (55) can then be found with a finite-
dimensional optimisation over the $t_{i}$ (see section 5.4 for implementation
details). If $u^{+}:L^{\prime}\rightarrow W^{1,p}_{0}(M)$ is a _peak
selection_ of $J$, meaning

$u^{+}(v)\in\arg\max\\{J(u)\mid u\in[L,v]_{S}\\},$

it is thus sufficient to find

$\displaystyle v_{2}:=\arg\min\\{J(u^{+}(v))\mid v\in L^{\prime}\\}$ (57)

to get a second eigenfunction $u_{p}^{(2)}:=u^{+}(v_{2})$, i.e.
$J(u_{p}^{(2)})=\lambda_{p}^{(2)}$. The minimisation problem (57) can be
solved with a similar descent scheme as Algorithm 1: The crucial observation
is that $-d$ can be used as a descent direction for $v\mapsto J(u^{+}(v))$
[36, Lemma 2.5]. Note that this is subtly different from the obvious
application of (48) which would be that
${\frac{d}{dt}J(u^{+}(v)-td)|_{t=0}<0}$, as opposed to
${\frac{d}{dt}J(u^{+}(v-td))|_{t=0}<0}$. The corresponding scheme is given in
Algorithm 2.

Algorithm 2 Computation of a second eigenfunction

  1. 1.

Compute the first eigenfunction $u^{(1)}_{p}$ with Algorithm 1.

  2. 2.

Set $k=1$, choose an initial guess $v^{[1]}\in W^{1,p}_{0}(M)$ and a parameter
$\alpha>0$.

  3. 3.

While $\big{\|}\operatorname{grad}J(u^{+}(v^{[k]}))\big{\|}_{p}$ is too large
do

    1. (a)

Compute a descent direction $d^{[k]}:=-{\rm
sign}(t_{2}(v^{[k]}))\operatorname{grad}J(u^{+}(v^{[k]}))$.

    2. (b)

For $s\in\mathbb{R}$ define

$v^{[k+1]}(s):=\frac{v^{[k]}+sd^{[k]}}{\|v^{[k]}+sd^{[k]}\|_{p}}.$

    3. (c)

Set $s^{*}:=\alpha/\max(1,\|d^{[k]}\|_{p})$ and check whether

$J(u^{+}(v^{[k+1]}(s^{*})))\leq
J(u^{+}(v^{[k]}))-\frac{s^{*}}{4}|t_{2}(v^{[k]})|\|\nabla d^{[k]}\|_{2}^{2}$

is fulfilled. If it is not, keep halving $s^{*}$ until it is.

    4. (d)

Set $v^{[k+1]}:=v^{[k+1]}(s^{*})$.

    5. (e)

Increment $k\to k+1$.

    6. (f)

If $s^{*}$ is smaller than some tolerance, go to Step 4.

  4. 4.

Return $u_{p}^{(2)}:=u^{+}(v^{[k]})$.

##  5 Experiments

We perform a series of experiments to illustrate how the eigenfunctions of the
dynamic $p$-Laplacian – and in particular their level sets with optimal
Cheeger ratio – compare to those of the $2$-Laplacian. We choose a variety of
domains: fully periodic, partially periodic, and non-periodic, and a variety
of dynamics with differing levels of symmetry. We find that the Cheeger ratios
of many superlevel sets of the first eigenfunction improve considerably for
$p\rightarrow 1$, indicating that the dynamic $p$-Laplacian behaves similarly
to 2.4(b). We also find that the level set of the leading eigenfunction of the
$p$-Laplacian with _least Cheeger ratio_ is relatively insensitive to $p$ for
$1.3\leq p\leq 2$; this finding is replicated for the dynamic $p$-Laplacian.
The code for replicating our results can be found in the accompanying package
DynamicPLaplacian.jl.

###  5.1 The static unit square

We first look at $M=[0,1]^{2}$ and $T=\mathrm{id}$. We get the functional
corresponding to the classical Cheeger problem

$J(u)=\frac{\|\nabla u\|_{p}^{p}}{\|u\|_{p}^{p}}.$

The theory predicts that for $p\rightarrow 1$ the first eigenvalue
$\lambda_{p}$ converges to the Cheeger constant

$h(M)=\inf_{A\subset M}\frac{\ell_{1}(\partial
A)}{\ell_{2}(A)}=\frac{4-\pi}{2-\sqrt{\pi}}\approx 3.772,$

(see section 2.2), while the Cheeger set is a square with rounded corners of
radius $1/h(M)\approx 0.265$.

We approximate the first eigenfunction $u_{p}$ of $-\Delta_{p}$ with Dirichlet
boundary conditions and investigate how the Cheeger ratio of its level sets
changes with $p$. For the discretisation of the eigenproblem (21) we use a
finite element method with linear triangular Lagrange elements on a $100\times
100$ grid, implemented using the package Gridap.jl [5]. For the stopping
criterion in Algorithm 2, we use $\|{\operatorname{grad}}^{D}\
{J}^{D}\|_{p}\leq 10^{-3}$. Note that in this example there is no dynamics,
i.e. we set $T:=id$, and thus
${\operatorname{grad}}^{D}{J}^{D}=\operatorname{grad}J$. While for $p$ close
to $2$ this is generally achieved in a few iterations, the number of
iterations grows quickly as $p$ gets close to $1$; see Figure 4. We restrict
the presentation of results to the range of $p$ where the convergence
criterion was achieved, namely $1.3\leq p\leq 2.0$.

Figure 4:  Number of iterations needed to achieve the convergence criterion
${\|\operatorname{grad}J\|_{p}\leq 10^{-3}}$ (or
${\|{\operatorname{grad}}^{D}{J}^{D}\|_{p}\leq 10^{-3}}$ respectively). For
details see sections 5.1 (static unit square), 5.2 (transitory double gyre),
5.4 (standard map) and 5.3 (cylinder flow).

In Figure 5, we plot the leading eigenfunction of $\Delta_{p}$ for three
values of $p$ (upper row) and their level sets (lower row). The level sets
were obtained by the marching squares algorithm implemented in Contour.jl
[11]. The first eigenfunction $u_{p}$ does indeed get “flatter” with smaller
$p$, and the level sets get closer to each other.

|  |
---|---|---
|  |
Figure 5: The first eigenfunction of $\Delta_{p}$ for $p=2.0$ (left), $p=1.6$
(centre) and $p=1.3$ (right) on the unit square.

We next compute the Cheeger ratio for specific level sets of the first
eigenfunction. For the computation of the area of the superlevel sets we used
the well-known formula sometimes called the _shoelace formula_ (see e.g. [32,
3.5 (6.)]. If the approximation of some eigenfunction is negative somewhere,
it can happen that a level set corresponding to a negative value has a very
small Cheeger ratio. We thus only consider level sets corresponding to
positive values in such cases. We observe in Figure 6 that the Cheeger ratio
of every single level set improves with decreasing $p$, while the Cheeger
ratios of the best level sets barely change with varying $p$. To quantify the
improvement, we investigate the median of the Cheeger ratio of the level sets
with respect to the uniform distribution on the range of the eigenfunction.

|  |
---|---|---
Figure 6: Cheeger ratio of level sets of the first eigenfunction of
$\Delta_{p}$ (left, closeup in the middle) on the unit square and statistics
of the ratio when choosing the level set randomly (right). Note that the ratio
does not depend monotonically on the level.

Finally, we check the best level set for different $p$ and compare it to the
boundary of the actual Cheeger set in Figure 7.

---
Figure 7: The level set with the smallest Cheeger ratio for different $p$
versus the actual Cheeger set (left) on the unit square. Closeup of the lower
left corner of the square (right).

###  5.2 The transitory double gyre

As a first example with dynamics we use the transitory double gyre introduced
in [27]. We consider the time one flow map of the non-autonomous differential
equation $(\dot{x},\dot{y})=(\partial_{y}\psi,-\partial_{x}\psi)$ on
$M=[0,1]^{2}$ defined by the stream function

$\displaystyle\psi(x,y,t)=(1-s(t))\psi_{P}(x,y)+s(t)\psi_{F}(x,y)$

with $\psi_{P}(x,y)=\sin(2\pi x)\sin(\pi y)$ and $\psi_{F}(x,y)=\sin(\pi
x)\sin(2\pi y)$ and

$s(t)=\begin{cases}0&\text{ for }t<0,\\\ t^{2}(3-2t)&\text{ for }t\in[0,1],\\\
1&\text{ for }t>1.\end{cases}$

In this flow, there are two vertically elongated gyres next to each other at
$t=0$ and these smoothly transition via rotation by $90$ degrees to two
horizontally elongated ones at ${t=1}$. We integrate the flow with the
Tsitouras 5/4 Runge–Kutta method implemented in the package
DifferentialEquations.jl [30] using a tolerance of $10^{-7}$ and approximate
the Cauchy–Green tensor with automatic differentiation, for which we use
ForwardDiff.jl [31].

The dynamic $p$-Laplacian, the first eigenfunction ${u_{p}}^{D}$, and its
level sets are approximated as in the previous example. Note that the dynamic
$p$-Laplacian acts on functions at the initial time $t=0$. As in the static
case, $u_{p}$ becomes flatter at its extremum for decreasing $p$; see Figure
8. In Figure 9, individual level sets and their statistics are shown, as well
as the optimal level sets for three different values of $p$. The overall shape
of the optimal set does not vary significantly with $p$.

|  |
---|---|---
|  |
Figure 8: The first eigenfunction of ${\Delta}^{D}_{p}$ for the transitory
double gyre and $p=2.0$ (left), $p=1.6$ (centre) and $p=1.3$ (right).

Figure 9: The dynamic Cheeger ratio of level sets of the first eigenfunction
of ${\Delta_{p}}^{D}$ for the transitory double gyre (top left), statistics of
the ratio when choosing the level set randomly (top right) and the level set
with the lowest Cheeger ratio for different $p$ (bottom left, closeup bottom
right).

###  5.3 The cylinder flow

As a more complex example we look at the flow of the vector field

$\displaystyle\dot{x}(t)$ $\displaystyle=c-A(t)\sin(x-\nu
t)\cos(y)+\varepsilon\Gamma(g(x,y,t))\sin(t/2)$ $\displaystyle\dot{y}(t)$
$\displaystyle=A(t)\cos(x-\nu t)\sin(y)$

on the domain $M=2\pi S^{1}\times[0,\pi]$, where $A(t)=1+\sin(2\sqrt{5}t)/8$,
$\Gamma(\psi)=1/(\psi^{2}+1)^{2}$, $g(x,y,t)=\sin(x-\nu t)\sin(y)+y/2-\pi/4$,
$c=0.5$, $\nu=0.5$, and $\varepsilon=0.25$, cf. [15, Section B].

As the flow time we choose $T=40$; the time integration is performed as in the
transitory double gyre example. The dynamic $p$-Laplacian is discretised on a
$200\times 100$ grid. Its leading eigenfunction and its level sets are
constructed as in the previous examples (see section 5.1 for details). As the
domain is periodic in one coordinate, the shoelace formula for the area of a
polygon does not hold anymore, and we instead determine the area of a
superlevel set by counting the number of vertices of a regular $200\times 100$
grid that lie in it.

The leading eigenfunctions $u^{(1)}_{p}$ of ${\Delta_{p}}^{D}$ for
$p\in\\{2.0,1.6,1.3\\}$ are depicted in Figure 10, and the dynamic Cheeger
ratios of the superlevel sets of the $u^{(1)}_{p}$ and their statistics are
shown in Figure 11. For the respective superlevel sets with minimal dynamic
Cheeger ratio see Figure 12.

|  |
---|---|---
|  |
Figure 10:  First eigenfunction of ${\Delta}^{D}_{p}$ for the cylinder flow
and $p=2.0$ (left), $p=1.6$ (centre) and $p=1.3$ (right).

Figure 11:  Dynamic Cheeger ratio of level sets of the first eigenfunction
(left) and statistics of the ratio when choosing the level set randomly
(right) for the cylinder flow.

Figure 12:  Level set with the smallest Cheeger ratio for different $p$ for
the cylinder flow (on the right a closeup of the region indicated in the left
plot).

###  5.4 The standard map

As an exploratory example, we consider the standard map on the flat $2$-torus
$\mathbb{T}^{2}$

$T(x,y)=(x+y+a\sin(x),y+a\sin(x))\;\;(\mathrm{mod}\;2\pi)$

with $a=0.971635$ as in [15]. Since the boundary is empty here, the first
eigenfunction $u_{p}^{(1)}$ is constant. We thus consider the second
eigenfunction $u_{p}^{(2)}$ and use Algorithm 2 to compute it. In order to
ensure solvability of eq. 54 we add a zero-mean condition $\int_{M}v=0$ in
step 3a. For convenience, we also add a zero-mean condition
$\int_{M}{d}^{D}=0$ in order to ensure uniqueness of the solution ${d}^{D}$.
For boundaryless domains, the relevant dynamic Cheeger constant is not (6),
but instead the quantity

${h}^{D}_{\rm Neum}(M,T):=\inf_{A\subset M}\frac{\ell_{d-1}(\partial
A)+\ell_{d-1}(\partial TA)}{2\min\\{\ell_{d}(A),\ell_{d}(M\backslash A)\\}},$

see [14]. We compute areas as in the previous example. The low-dimensional
optimisation problem that arises in the computation of $u^{+}(v)$ in Algorithm
2 was solved to local optimality with the L-BFGS method constrained to the
unit sphere in $\mathbb{R}^{2}$, as described in [1] and implemented in
Optim.jl [26]. As recommended in [36], we use the $t_{i}$ from the last
iteration to initialise the optimisation for the next iteration to promote
continuity of $u^{+}$. Otherwise, the computations have been done as in
section 5.1 and section 5.2. As vertical translations and $T$ are volume-
preserving and their Jacobian does not depend on $y$, ${F}^{D},G$ and thus
${J}^{D}$ are invariant under vertical translations. To deal with the
resulting non-uniqueness of eigenfunctions, we shift an eigenfunction such
that its maximum has $y$-coordinate $\pi$; this makes the eigenfunctions more
easily comparable.

In Figure 13, the second eigenfunction is shown for different $p$. We again
observe that the eigenfunction becomes flatter around its maximum as $p$
decreases. Similarly to the previous examples, the statistics of the Cheeger
ratio ${h}^{D}_{\rm Neum}(M,T)$ are given in Figure 14.

|  |
---|---|---
|  |

Figure 13: Second eigenfunction of ${\Delta}^{D}_{p}$ for the standard map and
$p=2.0$ (left), $p=1.6$ (centre) and $p=1.3$ (right)

.

Figure 14: Dynamic Cheeger ratio of level sets of the second eigenfunction of
${\Delta_{p}}^{D}$ for the standard map (left), statistics of the ratio when
choosing the level set randomly (center) and the level set with the lowest
Cheeger ratio for different $p$ (right).

##  6 Conclusions and outlook

We introduced a dynamic version of the nonlinear $p$-Laplace operator –
denoted ${\Delta}^{D}_{p}$ – suitable for the spectral analysis of dynamical
systems. In 3.1 we proved the existence of a positive leading eigenvalue
${\lambda}^{D}_{p}$ for $-{\Delta}^{D}_{p}$ with homogeneous Dirichlet
boundary conditions. Associated to the dynamic spectral geometry of the
dynamical system is the dynamic Cheeger constant ${h}^{D}(M,T)$ defined in
(6), which quantifies how well the boundary of a carefully chosen subset in
$M$ can resist growth under evolution of the dynamics. In Equation 22 we
developed a Cheeger-type inequality, relating ${\lambda}^{D}_{p}$ and the
dynamic Cheeger constant for $p\geq 1$. 3.5 showed that this inequality
becomes increasingly tight as $p$ approaches 1, and in fact that $\lim_{p\to
1}\lambda^{D}_{p}=h^{D}(M,T)$.

We then turned our attention to the numerical approximation of the dynamic
$p$-Laplacian and its eigenfunctions. We proposed an algorithm – modifying a
similar algorithm for the classical $p$-Laplacian due to Yao & Zhou [36] – to
solve the nonlinear eigenproblems for $1<p<2$, and compared the algorithm
performance and solution quality for $p$ between 1 and 2 on a variety of
examples. We found that (i) the level sets of the leading eigenfunctions
tended to increasingly concentrate around the boundary of the optimal coherent
set as $p\to 1$ and (ii) the optimal level set for the dynamic 2-Laplacian was
often not very different from the optimal level set of the dynamic
$p$-Laplacian for $p$ close to 1. This first finding suggests that the
“sharper” eigenfunctions of the dynamic $p$-Laplacian may be easier to use to
find coherent sets than those of the dynamic $2$-Laplacian, including using
multiple eigenfunctions post-processed via feature-separation algorithms such
as [17]. The second finding suggests that not too much is lost in terms of
solution quality (where coherent set identification is quantified via
$h^{D}(M,T)$) by using eigenfunctions of the dynamic 2-Laplacian.

There are numerous aspects of this work that could be addressed in the future.
One practical motivation for the current work was that (first) eigenfunctions
of $\Delta_{p}$ converge to characteristic functions of the Cheeger set.
Although our results deliver positive empirical evidence, such a result for
the dynamic $p$-Laplacian ${\Delta}^{D}_{p}$ has yet to be proven. Also, the
uniqueness and positivity of the first eigenfunction of ${\Delta}^{D}_{p}$ is
not yet established, as discussed in 3.1.

Finally, we have considered homogeneous Dirichlet boundary conditions, which
is consistent with looking for sets $A\subset M$ in the isoperimetric problem
(1) that have a compact closure in $M$. Neumann boundary conditions, on the
other hand, translate to problems where $\partial A$ can touch $\partial M$.
Results concerning Neumann boundary conditions for $\Delta_{p}$ are less well
developed, but extensions of our results to the dynamic $p$-Laplacian under
Neumann boundary conditions are natural to pursue.

## Acknowledgements

AD gratefully acknowledges support through the Bavarian Ministry of Science
and Art. The research of GF is partially supported by an Australian Research
Council Discovery Project, a Universities Australia Australia-Germany Joint
Research Cooperation Scheme, and an Einstein Visiting Fellowship to the Freie
Universität Berlin. GF is grateful to the Departments of Mathematics at the
Freie Universität Berlin and the University of Bayreuth for their generous
hospitality. PK has been partially supported by the Deutsche
Forschungsgemeinschaft (DFG) through grant CRC 1114 “Scaling Cascades in
Complex Systems”, Project Number 235221301, Project A08 “Characterization and
prediction of quasi-stationary atmospheric states”.

##  Appendix A Appendix

###  A.1 On the min-max theorem

######  Theorem A.1 (Courant–Fischer).

Let $A\in\mathbb{R}^{d\times d}$ be symmetric. The eigenvalues
$\lambda^{(1)}\leq\lambda^{(2)}\leq\cdots\leq\lambda^{(d)}$ of $A$ are given
by

$\lambda^{(k)}=\min\left\\{\max\left\\{\langle u,Au\rangle\mid u\in V\cap
S\right\\}\mid V\in V_{k}\right\\},$

where $V_{k}$ is the set of $k$-dimensional subspaces of $\mathbb{R}^{d}$ and
$S$ is the unit sphere in $\mathbb{R}^{d}$.

######  Proposition A.1.

In the setting of Theorem A.1, if $u_{k}$ denotes the
$\|\cdot\|_{2}$-normalised eigenvector at $\lambda^{(k)}$,
$L_{k}=\textrm{span}\left(u_{1},\ldots,u_{k}\right)=:[u_{1},\ldots,u_{k}]$,
and $L_{k}^{\prime}$ some fixed complement of $L_{k}$, then

$\lambda^{(k)}=\min\left\\{\max\left\\{\langle u,Au\rangle\mid
u\in[L_{k-1},v]\cap S\right\\}\mid v\in L_{k-1}^{\prime}\right\\}.$

Moreover, if $v_{k}$ denotes a minimiser of the above expression, then

$u_{k}\in\arg\max\big{\\{}\langle u,Au\rangle\mid u\in[L_{k-1},v_{k}]\cap
S\big{\\}}.$

###### Proof.

In the following computation, line by line, we make use of the following
facts:

  1. 1.

The Courant–Fischer theorem.

  2. 2.

For $v\in L_{k-1}^{\prime}$ one has $[L_{k-1},v]\in V_{k}$.

  3. 3.

Estimating the minimum from above by setting $v=v_{k}\in L_{k-1}^{\prime}$
such that $[L_{k-1},v_{k}]=L_{k}$.

  4. 4.

Using that the eigenvectors $u_{i}$ are mutually orthogonal and setting
$u=\sum_{i=1}^{k}t_{i}u_{i}$ with $\sum_{i=1}^{k}t_{i}^{2}=1$.

The computation

$\displaystyle\lambda^{(k)}$ $\displaystyle=\min\left\\{\max\left\\{\langle
u,Au\rangle\mid u\in V\cap S\right\\}\mid V\in V_{k}\right\\}$
$\displaystyle\leq\min\left\\{\max\left\\{\langle u,Au\rangle\mid
u\in[L_{k-1},v]\cap S\right\\}\mid v\in L_{k-1}^{\prime}\right\\}$
$\displaystyle\leq\max\left\\{\langle u,Au\rangle\mid u\in L_{k}\cap
S\right\\}$
$\displaystyle=\max\left\\{\textstyle\sum_{i=1}^{k}\lambda^{(i)}t_{i}^{2}\,\mid\,\sum_{i=1}^{k}t_{i}^{2}=1\right\\}\leq\left\\{\textstyle\sum_{i=1}^{k}\lambda^{(k)}t_{i}^{2}\,\mid\,\sum_{i=1}^{k}t_{i}^{2}=1\right\\}=\lambda^{(k)}$

yields the claim. Equality holds in the last inequality for $t_{k}=1$ and
$t_{i}=0$ for $i\neq k$, showing that a maximiser is indeed $u=u_{k}$. ∎

###  A.2 A lemma for the proof of Equation 22

######  Lemma A.1.

Let $M\in\mathbb{R}^{d}$ compact and $w\in C^{\infty}_{0}(M)$ and define

$\displaystyle A(t):=\\{x\mid w(x)>t\\}$ (58) $\displaystyle
B(t):=\\{x\mid|w(x)|>t\\}.$ (59)

  1. a)

For almost all $t$ we have $\partial A(t)=w^{-1}(t)$

  2. b)

For almost all $t$ we have $\partial B(t)=\partial A(t)\,\dot{\cup}\,\partial
A(-t)$

  3. c)

The following integrals coincide:

$\int_{-\infty}^{\infty}\ell_{d-1}(\partial A(t))\
dt=\int_{0}^{\infty}\ell_{d-1}(\partial B(t))dt.$

In all statements the boundary is taken within $M$.

###### Proof.

The inclusion $\partial A(t)\subseteq w^{-1}(t)$ follows for all $t$ by
continuity of $w$. On the other hand let $x\in w^{-1}(t)$. We may assume
$\nabla w(x)\neq 0$ as by Sard’s theorem $w^{-1}(t)$ contains no critical
points for almost all $t$. We may also assume $t\neq 0$ (as $\\{0\\}$ is a
null set), implying $x\in M^{\circ}$ and in particular $x+s\nabla w(x)\in M$
for small $s$. Now

$w(x+s\nabla w(x))-w(x)=s|\nabla w(x)|^{2}+\mathcal{O}(s^{2})$ (60)

is positive for small $s>0$ and negative for small $s<0$ and it follows that
$x\in\partial A(t)$. Hence, $\partial A(t)\supseteq w^{-1}(t)$ is fulfilled,
too and we have shown (a). Now this implies that for almost all $t$ we have

$\partial B(t)=\partial A(t)\,\dot{\cup}\,\partial A(-t)$ (61)

in the following way: if $x\in\partial B(t)$ then, again by continuity of $w$
we must have $w(x)=\pm t$, which for almost all $t$ implies $x\in\partial
A(t)\cup\partial A(-t)$ by (a). Conversely, if $x\in\partial A(\pm t)$ then
almost surely $x\in\partial B(t)$, as there are values bigger and smaller than
$w(x)$ attained arbitrarily close to $x$ by the same construction as in (60).
Finally, the sets $w^{-1}(t)$ and $w^{-1}(-t)$ are obviously disjoint for
$t\neq 0$ and so are the sets $\partial A(t)$ and $\partial A(-t)$ for almost
all $t$. This finishes the proof of (b). Finally, we can use this to conclude

$\displaystyle\int_{0}^{\infty}\ell_{d-1}(\partial B(t))dt$
$\displaystyle=\int_{0}^{\infty}\ell_{d-1}(\partial A(t))+\ell_{d-1}(\partial
A(-t))\ dt$ (62) $\displaystyle=\int_{0}^{\infty}\ell_{d-1}(\partial
A(t))+\int_{-\infty}^{0}\ell_{d-1}(\partial A(t))\ dt$ (63)
$\displaystyle=\int_{-\infty}^{\infty}\ell_{d-1}(\partial A(t)).$ (64)

∎

## References

  * [1] P.-A. Absil, R. Mahony, and R. Sepulchre.  Optimization algorithms on matrix manifolds.  In Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2009.
  * [2] R. A. Adams and J. J. Fournier.  Sobolev spaces.  Elsevier, 2003.
  * [3] H. Aref.  Stirring by chaotic advection.  Journal of Fluid Mechanics, 143:1–21, 1984.
  * [4] H. Aref, J. R. Blake, M. Budišić, S. S. Cardoso, J. H. Cartwright, H. J. Clercx, K. El Omari, U. Feudel, R. Golestanian, E. Gouillart, G. F. van Heijst, T. S. Krasnopolskaya, Y. L. Guer, R. S. MacKay, V. V. Meleshko, G. Metcalfe, I. Mezić, A. P. S. de Moura, O. Piro, M. F. M. Speetjens, R. Sturman, J.-L. Thiffeault, and I. Tuval.  Frontiers of chaotic advection.  Reviews of Modern Physics, 89(2):025007, 2017.
  * [5] S. Badia and F. Verdugo.  Gridap: An extensible finite element toolbox in Julia.  Journal of Open Source Software, 5(52):2520, 2020.
  * [6] M. Badiale and E. Serra.  Semilinear Elliptic Equations for Beginners.  Springer, London, 2011.
  * [7] R. Banisch and P. Koltai.  Understanding the geometry of transport: Diffusion maps for Lagrangian trajectory data unravel coherent sets.  Chaos, 27(3):035804, 2017.
  * [8] I. Chavel.  Eigenvalues in Riemannian geometry.  Academic Press, 1984.
  * [9] I. Chavel.  Isoperimetric Inequalities: Differential Geometric and Analytic Perspectives.  Number 145 in Cambridge tracts in mathematics. Cambridge University Press, Cambridge ; New York, 2001.
  * [10] J. Cheeger.  A Lower Bound for the Smallest Eigenvalue of the Laplacian.  In R. C. Gunning, editor, Problems in Analysis: A Symposium in Honor of Salomon Bochner (PMS-31), pages 195–200. Princeton University Press, 1970.
  * [11] D. Darakananda and T. Lycke.  Contour.jl.  https://github.com/JuliaGeometry/Contour.jl.  Accessed: June 2, 2023.
  * [12] H. Federer and W. H. Fleming.  Normal and integral currents.  Annals of Mathematics, pages 458–520, 1960.
  * [13] X. Feng and A. Prohl.  Analysis of total variation flow and its finite element approximations.  ESAIM: Mathematical Modelling and Numerical Analysis, 37(3):533–556, 2003.
  * [14] G. Froyland.  Dynamic isoperimetry and the geometry of Lagrangian coherent structures.  Nonlinearity, 28(10):3587–3622, 2015.
  * [15] G. Froyland and O. Junge.  On fast computation of finite-time coherent sets using radial basis functions.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(8):087409, 2015.
  * [16] G. Froyland and O. Junge.  Robust FEM-based extraction of finite-time coherent sets using scattered, sparse, and incomplete trajectories.  SIAM Journal on Applied Dynamical Systems, 17(2):1891–1924, 2018\.
  * [17] G. Froyland, C. P. Rock, and K. Sakellariou.  Sparse eigenbasis approximation: Multiple feature extraction across spatiotemporal scales with application to coherent set identification.  Communications in Nonlinear Science and Numerical Simulation, 77:81–107, 2019.
  * [18] G. Haller and F. J. Beron-Vera.  Coherent Lagrangian vortices: The black holes of turbulence.  Journal of Fluid Mechanics, 731:R4, 2013.
  * [19] J. Horák.  Numerical investigation of the smallest eigenvalues of the $p$-Laplace operator on planar domains.  Electronic Journal of Differential Equations, 2011(132):1–30, 2011\.
  * [20] B. Kawohl and V. Fridman.  Isoperimetric estimates for the first eigenvalue of the $p$-Laplace operator and the Cheeger constant.  Commentationes Mathematicae Universitatis Carolinae, 44:659–667, 2003.
  * [21] L. Lefton and D. Wei.  Numerical approximation of the first eigenpair of the $p$-Laplacian using finite elements and the penalty method.  Numerical Functional Analysis and Optimization, 18(3-4):389–399, 1997.
  * [22] G. P. Leonardi.  An overview on the Cheeger problem.  In New trends in shape optimization, volume 166, pages 117–139. Birkhäuser, Cham, 2015.
  * [23] P. Lindqvist.  On the equation $\operatorname{div}(|\nabla u|^{p-2}\nabla u)+\lambda|u|^{p-2}u=0$.  Proceedings of the American Mathematical Society, pages 157–164, 1990.
  * [24] P. Lindqvist.  A nonlinear eigenvalue problem.  Topics in mathematical analysis, 3:175–203, 2008.
  * [25] V. G. Maz’ya.  Classes of domains and imbedding theorems for function spaces.  In Doklady Akademii Nauk, volume 133, pages 527–530. Russian Academy of Sciences, 1960.
  * [26] P. K. Mogensen and A. N. Riseth.  Optim: A mathematical optimization package for Julia.  Journal of Open Source Software, 3(24):615, 2018.
  * [27] B. A. Mosovsky and J. D. Meiss.  Transport in transitory dynamical systems.  SIAM Journal on Applied Dynamical Systems, 10(1):35–65, 2011.
  * [28] E. Parini.  The second eigenvalue of the $p$-Laplacian as $p$ goes to 1.  International Journal of Differential Equations, 2010:1–23, 2010\.
  * [29] E. Parini.  An introduction to the Cheeger Problem.  Surv. Math. Appl., 6:9–21, 2011.
  * [30] C. Rackauckas and Q. Nie.  Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia.  Journal of Open Research Software, 5(1), 2017.
  * [31] J. Revels, M. Lubin, and T. Papamarkou.  Forward-mode automatic differentiation in Julia.  arXiv:1607.07892 [cs.MS], 2016.
  * [32] L. Råde and B. Westergren.  Mathematics Handbook for Science and Engineering.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.
  * [33] C. G. Simader.  On Dirichlet’s Boundary Value Problem.  Lecture Notes in Mathematics, Vol. 268, 1972.
  * [34] G. Strang and G. J. Fix.  An analysis of the finite element method.  Prentice-Hall Series in Automatic Computation. Englewood Cliffs, N.J.: Prentice-Hall, Inc. XIV, 306, 1973.
  * [35] R. Sturman, J. M. Ottino, and S. Wiggins.  The mathematical foundations of mixing: the linked twist map as a paradigm in applications: micro to macro, fluids to solids, volume 22.  Cambridge University Press, 2006.
  * [36] X. Yao and J. Zhou.  Numerical methods for computing nonlinear eigenpairs: Part I. Iso-homogeneous cases.  SIAM Journal on Scientific Computing, 29(4):1355–1374, 2007.