text:
\begin{document}

\title[Cancellable elements of the lattice of semigroup varieties]{Cancellable elements of the lattice\\
of semigroup varieties: varieties satisfying\\
a permutational identity of length 3}

\author{B.~M.~Vernikov}

\address{Ural Federal University, Institute of Natural Sciences and Mathematics, Lenina 51, 620083 Ekaterinburg, Russia}

\email{bvernikov@gmail.com}

\date{}

\thanks{The work was partially supported by the Ministry of Education and Science of the Russian Federation (project 1.6018.2017/8.9) and by the Russian Foundation for Basic Research (grant No.\,17-01-00551).}

\begin{abstract}
We completely determine all semigroup varieties satisfiyng a permutational identity of length 3 that are cancellable elements of the lattice of all semigroup varieties. Using this result, we provide a series of new examples of semigroup varieties that are modular but not cancellable elements of this lattice.
\end{abstract}

\keywords{Semigroup, variety, lattice of varieties, permutational identity, cancellable element of a lattice}

\subjclass{Primary 20M07, secondary 08B15}

\maketitle

\section{Introduction and summary}
\label{intr}

There are a number of articles devoted to an examination of special elements in the lattice $\mathbb{SEM}$ of all semigroup varieties (see surveys~\cite[Section~14]{Shevrin-Vernikov-Volkov-09} and~\cite{Vernikov-15}).One can recall definitions of two types of special elements that appear below. An element $x$ of a lattice $\langle L;\vee,\wedge\rangle$ is called
\begin{align*}
&\text{\emph{modular} if}\enskip&&(\forall y,z\in L)\quad \bigl(y\le z\longrightarrow(x\vee y)\wedge z=(x\wedge z)\vee y\bigr),\\
&\text{\emph{cancellable} if}\enskip&&(\forall y,z\in L)\quad (x\vee y=x\vee z\ \&\ x\wedge y=x\wedge z\longrightarrow y=z).
\end{align*}
It is evident that every cancellable element of a lattice is modular. A valuable information about modular and cancellable elements in abstract lattices can be found in~\cite{Seselja-Tepavcevic-01}, for instance.

Several results about modular elements of the lattice $\mathbb{SEM}$ were provided in the papers~\cite{Jezek-McKenzie-93,Shaprynskii-12,Vernikov-07}, while cancellable elements of $\mathbb{SEM}$ were considered in~\cite{Gusev-Skokov-Vernikov-18+,Skokov-Vernikov-18+}. In particular, commutative semigroup varieties that are modular elements in $\mathbb{SEM}$ are completely determined by the author  in~\cite[Theorem~3.1]{Vernikov-07}. Further, it is verified by S.Gusev, D.Skokov and the author in~\cite[Theorem~1.1]{Gusev-Skokov-Vernikov-18+} that the properties to be modular and cancellable elements of $\mathbb{SEM}$ are equivalent in the class of commutative semigroup varieties. The question is naturally arised, whether this equivalence true for arbitrary semigroup varieties. The negative answer to this question is given in~\cite{Skokov-Vernikov-18+}. Recall that an identity of the form
\begin{equation}
\label{permut id}
x_1x_2\cdots x_n\approx x_{1\pi}x_{2\pi}\cdots x_{n\pi}
\end{equation}
where $\pi$ is a non-trivial permutation on the set $\{1,2,\dots,n\}$ is called a \emph{permutational} identity. The number $n$ is called a \emph{length} of the identity~\eqref{permut id}. Clearly, the commutative law is a permutational identity of length~2. In~\cite[Theorem~1.1]{Skokov-Vernikov-18+} (see Proposition~\ref{mod-permut-3} below) modular elements of the lattice $\mathbb{SEM}$ are described within the class of varieties satisfying a permutational identity of length~3. After that, among varieties satisfying the conclusion of this statement, a variety is found that is not a cancellable element in $\mathbb{SEM}$~\cite[Proposition~1.2]{Skokov-Vernikov-18+}.

Here we continue examinationes started in~\cite{Gusev-Skokov-Vernikov-18+,Skokov-Vernikov-18+}. Namely, we completely classify varieties that satisfy a  permutational identity of length~3 and are cancellable elements in $\mathbb{SEM}$. Comparison of this result and~\cite[Theorem~1.1]{Skokov-Vernikov-18+} permits to specify a number of new examples of modular but not cancellable elements in $\mathbb{SEM}$.

A semigroup variety is called a \emph{nil-variety} if it consists of nilsemigroups. Semigroup words unlike letters are written in bold. Two sides of identities we connect by the symbol~$\approx$, while the symbol~$=$ stands for the equality relation on the free semigroup. As usual, we write the pair of identities $x\mathbf{u\approx u}x\approx\mathbf u$ where the letter $x$ does not occur in the word \textbf u in the short form $\mathbf u\approx 0$ and refer to the expression $\mathbf u\approx 0$ as to a single identity.  We denote by \textbf T the trivial semigroup variety and by \textbf{SL} the variety of all semilattices.

The main result of the article is the following

\begin{theorem}
\label{canc-permut-3}
A semigroup variety $\mathbf V$ satisfying a permutational identity of length $3$ is a cancellable element in the lattice $\mathbb{SEM}$ if and only if $\mathbf V=\mathbf{M\vee N}$ where $\mathbf M$ is one of the  varieties $\mathbf T$ or $\mathbf{SL}$, while the variety $\mathbf N$ satisfies the following identities:
\begin{align}
\label{xyz=yxz=xzy}xyz&\approx yxz\approx xzy,\\
\label{xxy=0}x^2y&\approx 0.
\end{align}
\end{theorem}

It is easy to see that if a variety satisfies the identities~\eqref{xyz=yxz=xzy} then it satisfies all permutational identities of length~3. Thus, Theorem~\ref{canc-permut-3} shows that if a cancellable element of $\mathbb{SEM}$ satisfies one permutational identity of length~3 then it satisfies all identities of such a form. As we will seen below, the analog of this claim is true for permutational identities of arbitrary length (see Proposition~\ref{alternative for cancellable}).

The article consists of three sections. Section~\ref{prel} contains auxiliary results. In Section~\ref{proof} we verify Theorem~\ref{canc-permut-3}.

\section{Preliminaries}
\label{prel}

First of all, we reproduce the main result of the article~\cite{Gusev-Skokov-Vernikov-18+}.

\begin{proposition}[\textup{\!\!\cite[Theorem~1.1]{Gusev-Skokov-Vernikov-18+}}]
\label{mod-permut-3}
A semigroup variety $\mathbf V$ satisfying a permutational identity of length $3$ is a modular element in the lattice $\mathbb{SEM}$ if and only if $\mathbf V=\mathbf{M\vee N}$ where $\mathbf M$ is one of the  varieties $\mathbf T$ or $\mathbf{SL}$, while the variety $\mathbf N$ satisfies one of the following identity systems: \textup{(i)}~$xyz\approx zyx$, $x^2y\approx 0$; \textup{(ii)}~$xyz\approx yzx$, $x^2y\approx 0$; \textup{(iii)}~$xyz\approx yxz$, $xyzt\approx xzty$, $xy^2\approx 0$; \textup{(iv)}~$xyz\approx xzy$, $xyzt\approx yzxt$, $x^2y\approx 0$.\qed
\end{proposition}

Comparison of Theorem~\ref{canc-permut-3} and Proposition~\ref{mod-permut-3} allows us to provide many examples of modular but not cancellable elements in the lattice $\mathbb{SEM}$. In particular, the varieties given by each of the identity systems (i)--(iv)  listed in Proposition~\ref{mod-permut-3} have this property.

\begin{lemma}
\label{modular non-cancellable}
\!\!\textup{(\!\!\cite[Lemma~2.3]{Gusev-Skokov-Vernikov-18+})} Let $x$ be a modular but not cancellable element of a lattice $L$ and let $y$ and $z$ be different elements of $L$ such that $x\vee y=x\vee z$ and $x\wedge y=x\wedge z$. Then there is an element $x'\in L$ such that $x'\le x$, $x'\vee y=x'\vee z$, $x'\wedge y=x'\wedge z$ and $y\vee z=x'\vee y$.\qed
\end{lemma}

The following lemma immediately follows from Lemmas~2.1 and~2.8 of the article~\cite{Gusev-Skokov-Vernikov-18+}.

\begin{lemma}
\label{join with SL}
A semigroup variety $\mathbf X$ is a cancellable element of the lattice $\mathbb{SEM}$ if and only if the variety $\mathbf{X\vee SL}$ has this property.\qed
\end{lemma}

We denote by $\mathbf{SEM}$ the variety of all semigroups. The following claim gives a strong necessary condition for a semigroup variety to be a modular element in the lattice $\mathbb{SEM}$.

\begin{proposition}
\label{nec}
If $\mathbf V$ is a modular  element of the lattice $\mathbb{SEM}$ then either $\mathbf V=\mathbf{SEM}$ or $\mathbf V=\mathbf{M\vee N}$ where $\mathbf M$ is one of the varieties $\mathbf T$ or $\mathbf{SL}$, while  $\mathbf N$ is a nil-variety.\qed
\end{proposition}

This proposition was proved (in slightly weaker form and in some other terminology) in~\cite[Proposition~1.6]{Jezek-McKenzie-93}. A deduction of Proposition~\ref{nec} from~\cite[Proposition~1.6]{Jezek-McKenzie-93} was given explicitly in~\cite[Proposition~2.1]{Vernikov-07}. A direct and transparent proof of Proposition~\ref{nec} not depending on a technique from~\cite{Jezek-McKenzie-93} is provided in~\cite{Shaprynskii-12}. Note that an essetially stronger necessary condition for a semigroup variety to be a modular element in the lattice $\mathbb{SEM}$ than Proposition~\ref{nec} is given by the author in~\cite[Theorem~2.5]{Vernikov-07} but this stronger results will not be used below.

The following claim can be easily deduced from~\cite[Lemma~1]{Sapir-Sukhanov-81}.

\begin{lemma}
\label{split}
If a nil-variety of semigroups $\mathbf N$ satisfies an identity of the form
\begin{equation}
\label{x_1x_2...x_n=v}
x_1x_2\cdots x_n\approx\mathbf v
\end{equation}
then either $\mathbf N$ satisfies the identity
\begin{equation}
\label{nilpotence}
x_1x_2\cdots x_n\approx 0
\end{equation}
or the identity~\eqref{x_1x_2...x_n=v} is permutational.\qed
\end{lemma}

Let $n$ be a natural number. We denote by $S_n$ the full permutation group on the set $\{1,2,\dots,n\}$. If $1\le i\le n$ then we denote by $\Stab_n(i)$ the set of all permutations $\pi\in S_n$ with $i\pi=i$. Obviously, $\Stab_n(i)$ is a subgroup in $S_n$. Moreover, it is well known that $\Stab_n(i)$ is a maximal proper subgroup in $S_n$. Let $T$ be the trivial group, $T_{ij}$ be the group generated by the transposition $(ij)$, $C_{ijk}$ and $C_{ijk\ell}$ be the groups generated by the cycles $(ijk)$ and $(ijk\ell)$ respectively, $P_{ij,k\ell}$ be the group generated by the disjoint transpositions $(ij)$ and $(k\ell)$, $A_n$ be the alternative subgroup of $S_n$ and $V_4$ be the Klein four-group. The subgroup lattice of the group $G$ is denoted by $\Sub(G)$. We need to know the structure of the lattices $\Sub(S_3)$ and $\Sub(S_4)$.  It is generally known and easy to check that the first of these two lattices has the form shown in Fig.~\ref{Sub(S_3)}. Direct routine calculations allow to verify that the lattice $\Sub(S_4)$ is as shown in Fig.~\ref{Sub(S_4)}.

\begin{figure}[tbh]
\unitlength=1mm
\special{em:linewidth 0.4pt}
\linethickness{0.4pt}
\begin{picture}(73,29)
\put(6,14){\circle*{1.33}}
\put(21,14){\circle*{1.33}}
\put(36,4){\circle*{1.33}}
\put(36,24){\circle*{1.33}}
\put(51,14){\circle*{1.33}}
\put(66,14){\circle*{1.33}}
\gasset{AHnb=0}
\drawline(36,4)(6,14)(36,24)(21,14)(36,4)(51,14)(36,24)(66,14)(36,4)
\put(5,14){\makebox(0,0)[rc]{$T_{12}$}}
\put(20,14){\makebox(0,0)[rc]{$T_{13}$}}
\put(52,14){\makebox(0,0)[lc]{$T_{23}$}}
\put(67,14){\makebox(0,0)[lc]{$C_{123}$}}
\put(36,27){\makebox(0,0)[cc]{$S_3$}}
\put(36,1){\makebox(0,0)[cc]{$T$}}
\end{picture}
\caption{The lattice $\Sub(S_3)$}
\label{Sub(S_3)}
\end{figure}

\begin{figure}[tbh]
\unitlength=0.9mm
\special{em:linewidth 0.4pt}
\linethickness{0.4pt}
\begin{picture}(135,70)
\put(73,5){\circle*{1.4}}
\put(13,20){\circle*{1.4}}
\put(23,20){\circle*{1.4}}
\put(33,20){\circle*{1.4}}
\put(43,20){\circle*{1.4}}
\put(53,20){\circle*{1.4}}
\put(63,20){\circle*{1.4}}
\put(73,20){\circle*{1.4}}
\put(83,20){\circle*{1.4}}
\put(93,20){\circle*{1.4}}
\put(103,20){\circle*{1.4}}
\put(113,20){\circle*{1.4}}
\put(123,20){\circle*{1.4}}
\put(133,20){\circle*{1.4}}
\put(3,35){\circle*{1.4}}
\put(13,35){\circle*{1.4}}
\put(23,35){\circle*{1.4}}
\put(33,35){\circle*{1.4}}
\put(43,35){\circle*{1.4}}
\put(53,35){\circle*{1.4}}
\put(63,35){\circle*{1.4}}
\put(73,35){\circle*{1.4}}
\put(83,35){\circle*{1.4}}
\put(93,35){\circle*{1.4}}
\put(103,35){\circle*{1.4}}
\put(73,50){\circle*{1.4}}
\put(83,50){\circle*{1.4}}
\put(93,50){\circle*{1.4}}
\put(113,50){\circle*{1.4}}
\put(73,65){\circle*{1.4}}
\gasset{AHnb=0}
\drawline(73,5)(13,20)(3,35)(23,20)(73,5)(33,20)(3,35)(73,65)(73,5)(83,20)(103,35)(93,20)(73,5)(103,20)(113,50)(93,20)(73,5)(123,20)(113,50)(133,20)(73,5)(43,20)(23,35)(73,65)(13,35)(13,20)(73,50)(103,35)(83,50)(73,65)(93,50)(103,35)(113,50)(73,65)(33,35)(33,20)(93,50)(93,20)(63,35)(63,20)(23,35)(123,20)
\drawline(23,35)(23,20)(83,50)(83,20)(53,35)(53,20)(73,5)(63,20)(13,35)(53,20)(33,35)(43,20)(43,35)(73,20)(103,35)
\drawline(33,20)(33,35)(133,20)
\drawline(13,35)(113,20)(73,5)
\drawline(113,20)(113,50)
\drawline(3,35)(103,20)
\footnotesize
\put(113,53){\makebox(0,0)[cc]{$A_4$}}
\put(128,17){\makebox(0,0)[lc]{$C_{123}$}}
\put(87,38){\makebox(0,0)[rc]{$C_{1234}$}}
\put(118,17){\makebox(0,0)[lc]{$C_{124}$}}
\put(97,38){\makebox(0,0)[rc]{$C_{1243}$}}
\put(77,38){\makebox(0,0)[rc]{$C_{1324}$}}
\put(108,17){\makebox(0,0)[lc]{$C_{134}$}}
\put(98,17){\makebox(0,0)[lc]{$C_{234}$}}
\put(76,17){\makebox(0,0)[rc]{$P_{12,34}$}}
\put(86,17){\makebox(0,0)[rc]{$P_{13,24}$}}
\put(96,17){\makebox(0,0)[rc]{$P_{14,23}$}}
\put(73,68){\makebox(0,0)[cc]{$S_4$}}
\put(6,43){\makebox(0,0)[cc]{$\Stab_4(1)$}}
\put(16,48){\makebox(0,0)[cc]{$\Stab_4(2)$}}
\put(31,53){\makebox(0,0)[cc]{$\Stab_4(3)$}}
\put(46,58){\makebox(0,0)[cc]{$\Stab_4(4)$}}
\put(73,2){\makebox(0,0)[cc]{$T$}}
\put(43,17){\makebox(0,0)[cc]{$T_{12}$}}
\put(53,17){\makebox(0,0)[cc]{$T_{13}$}}
\put(63,17){\makebox(0,0)[cc]{$T_{14}$}}
\put(33,17){\makebox(0,0)[cc]{$T_{23}$}}
\put(23,17){\makebox(0,0)[cc]{$T_{24}$}}
\put(13,17){\makebox(0,0)[cc]{$T_{34}$}}
\put(104,32){\makebox(0,0)[cc]{$V_4$}}
\gasset{AHnb=1}
\drawline(4,42)(3,36)
\drawline(14,47)(13,36)
\drawline(29,52)(23,36)
\drawline(44,57)(33,36)
\end{picture}
\caption{The lattice  $\Sub(S_4)$}
\label{Sub(S_4)}
\end{figure}

We need the following

\begin{lemma}
\label{subgroups cancellable}
A subgroup $G$ of the group $S_n$ is a cancellable element of the lattice $\Sub(S_n)$ if and only if either $G=T$ or $G=S_n$.
\end{lemma}

\begin{proof}
If $n\le 2$ then $S_n$ does not contain subgroups differ from $T$ and $S_n$. If $n=3$ then the desirable conclusion immediately follows from Fig.~\ref{Sub(S_3)}. Let now $n\ge 4$ and $G$ be a non-singleton proper subgroup of $S_n$. Suppose that $G$ is a cancellable and therefore, modular element of $\Sub(S_n)$. If $n=4$ then $G\supseteq V_4$ by~\cite[Proposition~3.8]{Jezek-81}, while if $n\ge 5$ then $G\supseteq A_n$ by~\cite[Proposition~3.1]{Jezek-81}. Clearly, it suffices to verify that there are at least two complements to $G$ in $\Sub(S_n)$. Suppose that $G\supseteq A_n$. Then $G=A_n$ because $A_n$ is a maximal proper subgroup in $S_n$. Then all subgroups of the form $T_{ij}$ are complements to $G$ in $\Sub(S_n)$. It remains to consider the case when $n=4$, $V_4\subseteq G\subset S_4$ and $G\ne A_4$. Fig.~\ref{Sub(S_4)} implies that either $G=V_4$ or $G=V_4\vee P_{ij,k\ell}$ for some disjoint transpositions $(ij)$ and $(k\ell)$. If $G=V_4$ then all subgroups of the form $\Stab_4(i)$ are complements to $G$ in $\Sub(S_4)$. Finally, if $G=V_4\vee P_{ij,k\ell}$ then subgroups $T_{ik}$ and $T_{j\ell}$ are complements to $G$ in $\Sub(S_4)$.
\end{proof}

\section{The proof of Theorem~\ref{canc-permut-3}}
\label{proof}

For a semigroup variety \textbf V, we denote by $\Perm_n(\mathbf V)$ the set of all permutations $\pi\in S_n$ such that \textbf V satisfies the identity~\eqref{permut id}. Clearly, $\Perm_n(\mathbf V)$ is a subgroup in $S_n$. We need the following

\begin{lemma}
\label{modular from V to Perm_n(V)}
\!\!\textup{(\!\!\cite[Corollary~4.3]{Vernikov-07})} If a semigroup variety $\mathbf V$ is a modular element of the lattice $\mathbb{SEM}$ and $n$ is a positive integer then the group $\Perm_n(\mathbf V)$ is a modular element of the lattice $\Sub(S_n)$.\qed
\end{lemma}

A word $\mathbf u$ is called \emph{linear} if every letter occurs in $\mathbf u$ at most once. If $\mathbf u$ is a word then we denote by $\ell(\mathbf u)$ the length of $\mathbf u$ and by $\con(\mathbf u)$ the set of letters occurring in $\mathbf u$.  For brevity, we will denote the identity~\eqref{permut id} by $p_n[\pi]$. If $\mathbf w$ is a word, $\con(\mathbf w)=\{x_1,x_2,\dots,x_n\}$ and $\xi\in S_n$ then we denote by $\xi[\mathbf w]$ the word obtained from $\mathbf w$ by the substitution $x_i\mapsto x_{i\xi}$ for all $i=1,2,\dots,n$.

\begin{proposition}
\label{cancellable from V to Perm_n(V)}
If a semigroup variety $\mathbf V$ is a cancellable element of the lattice $\mathbb{SEM}$ and $n$ is a positive integer then the group $\Perm_n(\mathbf V)$ is a cancellable element of the lattice $\Sub(S_n)$.
\end{proposition}

\begin{proof}
Clearly, $\mathbf V$ is a modular  element of the lattice $\mathbb{SEM}$. By Proposition~\ref{nec}, either $\mathbf V=\mathbf{SEM}$ or $\mathbf V=\mathbf{M\vee N}$ where $\mathbf M$ is one of the varieties $\mathbf T$ or $\mathbf{SL}$, while  $\mathbf N$ is a nil-variety. It is evident that $\Perm_n(\mathbf{SEM})=T$ and $\Perm_n(\mathbf{SL\vee N})=\Perm_n(\mathbf N)$ for any $n$. Since $T$ is a cancellable element of $\Sub(S_n)$, we can assume that $\mathbf V=\mathbf N$. In particular, $\mathbf V$ is a nil-variety.

Put $V=\Perm_n(\mathbf V)$ for brevity. Suppose that $V$ is not a cancellable element of the lattice $\Sub(S_n)$. Then there are subgroups $X_1$ and $X_2$ of the group $S_n$ such that $V\vee X_1=V\vee X_2$ and $V\wedge X_1=V\wedge X_2$ but $X_1\ne X_2$. For $i=1,2$, we denote by $\mathbf X_i$ the variety given by the identity $x_1x_2\cdots x_{n+1}\approx 0$, all identities of the form $\mathbf w\approx 0$ where $\mathbf w$ is a word of length $n$ depending on $<n$ letters and all identities of the form $p_n[\pi]$ where $\pi\in X_i$. It is clear that $\mathbf X_1\ne\mathbf X_2$.

One can verify that none of the varieties $\mathbf V$, $\mathbf X_1$ and $\mathbf X_2$ satisfies the identity~\eqref{nilpotence}. Indeed, suppose that this identity holds in $\mathbf V$. Then $\mathbf V$ satisfies all permutational identities of length $n$, whence $V=S_n$. But then
$$
X_1=V\wedge X_1=V\wedge X_2=X_2,
$$
contradicting the choice of the groups $X_1$ and $X_2$. Suppose now that the identity~\eqref{nilpotence} holds in $\mathbf X_1$. Then $X_1=S_n\supseteq X_2$. Lemma~\ref{modular from V to Perm_n(V)} implies that $V$ is a modular element of the lattice $\Perm_n(\mathbf V)$. Therefore,
\begin{align*}
X_1={}&(V\vee X_1)\wedge X_1&&\text{by the absorbtion law}\\
={}&(V\vee X_2)\wedge X_1&&\text{because }V\vee X_1=V\vee X_2\\
={}&(V\wedge X_1)\vee X_2&&\text{because }X_1\supseteq X_2\text{ and }V\text{ is modular in }\Perm_n(\mathbf V)\\
={}&(V\wedge X_2)\vee X_2&&\text{because }V\wedge X_1=V\wedge X_2\\
={}&X_2&&\text{by the absorbtion law},
\end{align*}
contradicting the choice of the groups $X_1$ and $X_2$ again. Analogous arguments show that the identity~\eqref{nilpotence} is not satisfied by $\mathbf X_2$.

Let now $\mathbf{u\approx v}$ be an arbitrary identity that holds in $\mathbf{V\vee X}_1$. We are going to verify that this identity holds in $\mathbf{V\vee X}_2$. Since $\mathbf{u\approx v}$ holds in $\mathbf V$, it suffices to verify that it holds in $\mathbf X_2$.  The identity  $\mathbf{u\approx v}$ holds in $\mathbf V$ and $\mathbf X_1$. If $\ell(\mathbf u),\ell(\mathbf v)\ge n+1$ then $\mathbf u\approx 0\approx \mathbf v$ in $\mathbf X_2$. Thus, we can assume without loss of generality that $\ell(\mathbf u)\le n$. On the other hand, the definition of the variety $\mathbf X_1$ and the fact that $\mathbf{u\approx v}$ holds in $\mathbf X_1$ imply that $\ell(\mathbf u),\ell(\mathbf v)\ge n$. In particular, $\ell(\mathbf u)=n$.

Suppose that the word $\mathbf u$ is linear. By Lemma~\ref{split}, either the identity $\mathbf{u\approx v}$ is permutational or the variety $\mathbf{V\vee X_1}$ satisfies the identity $\mathbf u\approx 0$. But we have proved above that the second case is impossible. Therefore, $\mathbf{u\approx v}$ is an identity of the form $p_n[\pi]$. Since it holds both in $\mathbf V$ and $\mathbf X_1$, we have $\pi\in V\wedge X_1=V\wedge X_2$. Hence $\pi\in X_2$, and therefore the identity $\mathbf{u\approx v}$ holds in $\mathbf X_2$.

It remains to consider the case when the word $\mathbf u$ is non-linear. Then it depends on $<n$ letters. Therefore, $\mathbf u\approx 0$ in the varieties $\mathbf X_1$ and $\mathbf X_2$. If $\ell(\mathbf v)>n$ or $\mathbf v$ is a word of length $n$ depending on $<n$ letters then $\mathbf v\approx 0$ in $\mathbf X_2$. Then $\mathbf{u\approx v}$ holds in $\mathbf X_2$. Finally, if $\ell(\mathbf v)=n$ and $\mathbf v$ depends on $n$ letters then the word $\mathbf v$ is linear and we can complete our considerations by the same arguments as in the previous paragraph.

Thus, if the identity $\mathbf{u\approx v}$ holds in $\mathbf{V\vee X}_1$ then it holds in $\mathbf{V\vee X}_2$ too. This means that $\mathbf{V\vee X}_2\subseteq\mathbf{V\vee X}_1$. The inverse inclusion can be verified analogously, whence $\mathbf{V\vee X}_1=\mathbf{V\vee X}_2$.

Let now $\mathbf{u\approx v}$ be an arbitrary identity that holds in the variety $\mathbf{V\wedge X}_1$. We aim to verify that it holds in $\mathbf{V\wedge X}_2$. Let the sequence of words
\begin{equation}
\label{deduction}
\mathbf u=\mathbf u_0,\mathbf u_1,\dots,\mathbf u_k=\mathbf v
\end{equation}
be the shortest deduction of the identity $\mathbf{u\approx v}$ from the identities of the varieties $\mathbf V$ and $\mathbf X_1$. This means that, for any $i=0,1,\dots,k-1$, the identity $\mathbf u_i\approx\mathbf u_{i+1}$ holds in one of the varieties $\mathbf V$ and $\mathbf X_1$ and the situation when one of the varieties $\mathbf V$ and $\mathbf X_1$ satisfies the identities $\mathbf u_i\approx\mathbf u_{i+1}\approx\mathbf u_{i+2}$ for some $0\le i\le k-2$  is impossible.

Suppose that there is an index $i$ such that $\mathbf u_i$ a linear word of length $n$. If $i>0$ then Lemma~\ref{split} implies that either $\mathbf u_{i-1}$ is a linear word and $\con(\mathbf u_{i-1})=\con(\mathbf u_i)$ or one of the varieties $\mathbf V$ or $\mathbf X_1$ satisfies the identity~\eqref{nilpotence}. Analogously, if $i<k$ then either $\mathbf u_{i+1}$ is a linear word and $\con(\mathbf u_i)=\con(\mathbf u_{i+1})$ or one of the varieties $\mathbf V$ and $\mathbf X_1$ satisfies the identity~\eqref{nilpotence}. But we have seen above that the last identity fails in both  the varieties $\mathbf V$ and $\mathbf X_1$. Therefore, the words adjacent to $\mathbf u_i$ in the sequence~\eqref{deduction} are linear words of length $n$ depending on the same letters as $\mathbf u_i$. By the trivial induction, this means that all the words $\mathbf u_0$, $\mathbf u_1$, \dots, $\mathbf u_k$ are linear words of length $n$ depending on the same letters. We may assume without loss of generality that $\con(\mathbf u)=\{x_1,x_2,\dots,x_n\}$. There are permutations $\pi_0,\pi_1,\dots,\pi_{k-1}\in S_n$ such that $\mathbf u_i=\pi_i[\mathbf u_{i-1}]$ for each $i=1,2,\dots,k$. Clearly, $\pi_i\in V$ [respectively $\pi_i\in X_1$] whenever the identity $\mathbf u_{i-1}\approx\mathbf u_i$ holds in the variety $\mathbf V$ [respectively $\mathbf X_1$]. Put $\pi=\pi_0\pi_1\cdots\pi_{k-1}$. Then $\pi\in V\vee X_1=V\vee X_2$. Therefore, there are permutations $\sigma_0,\sigma_1,\dots,\sigma_{m-1}\in S_n$ such that $\pi=\sigma_0\sigma_1\cdots\sigma_{m-1}$ and, for all $i=0,1,\dots,m-1$, the permutation $\sigma_i$ lies in either $V$ or $X_2$. Put $\mathbf v_0=\mathbf u$ and $\mathbf v_i=\sigma_i[\mathbf u_{i-1}]$ for each $i=1,2,\dots,m$. Obviously, $\mathbf v_m=\mathbf v$ and, for any $i=0,1,\dots,m-1$, the identity $\mathbf v_i\approx\mathbf v_{i+1}$ holds in one of the varieties $\mathbf V$ or $\mathbf X_2$. Therefore, the identity $\mathbf{u\approx v}$ holds in the variety $\mathbf{V\wedge X}_2$.

Suppose now that there are no linear words of length $n$ among the words $\mathbf u_0$, $\mathbf u_1$, \dots, $\mathbf u_k$. The definition of the variety $\mathbf X_1$ shows that if this variey satisfies an identity $\mathbf p=\mathbf q$ then $\ell(\mathbf p),\ell(\mathbf q)\ge n$. This means that:
\begin{itemize}
\item[$\bullet$] $\ell(\mathbf u_i)\ge n$ for all $i=1,2,\dots,k-1$;
\item[$\bullet$] either $\ell(\mathbf u_0)\ge n$ or the identity $\mathbf u_0\approx\mathbf u_1$ holds in $\mathbf V$;
\item[$\bullet$] either $\ell(\mathbf u_k)\ge n$ or the identity $\mathbf u_{k-1}\approx\mathbf u_k$ holds in $\mathbf V$.
\end{itemize}
In view of we say above, if $\ell(\mathbf u_i)=n$ for some $0\le i\le k$ then the word $\mathbf u_i$ is non-linear. Then the sequence~\eqref{deduction} is a deduction of the identity $\mathbf{u\approx v}$ from the identities of the varieties $\mathbf V$ and $\mathbf X_2$. Thus, we have again that  the identity $\mathbf{u\approx v}$ holds in the variety $\mathbf{V\wedge X}_2$.

We prove that if the identity $\mathbf{u\approx v}$ holds in $\mathbf{V\wedge X}_1$ then it holds in $\mathbf{V\wedge X}_2$. Therefore, $\mathbf{V\wedge X}_2\subseteq\mathbf{V\wedge X}_1$. The inverse inclusion can be verified analogously, whence $\mathbf{V\wedge X}_1=\mathbf{V\wedge X}_2$.

Thus, $\mathbf{V\vee X}_1=\mathbf{V\vee X}_2$ and $\mathbf{V\wedge X}_1=\mathbf{V\wedge X}_2$. Since $\mathbf V$ is cancellable in $\mathbb{SEM}$, this implies that $\mathbf X_1=\mathbf X_2$. But then $X_1=X_2$, contradicting the choice of the groups $X_1$ and $X_2$.
\end{proof}

Proposition~\ref{cancellable from V to Perm_n(V)} and Lemma~\ref{subgroups cancellable} immediately imply the following

\begin{proposition}
\label{alternative for cancellable}
Let $\mathbf V$ be a cancellable element of the lattice $\mathbb{SEM}$ and $n$ be a natural number. If $\mathbf V$ satisfies some permutational identity of length $n$ then it satisfies all such identities.\qed
\end{proposition}

Proposition~\ref{alternative for cancellable} and~\cite[Theorem~4.5]{Vernikov-07} imply the following

\begin{corollary}
\label{non-linear = 0}
Let $\mathbf V$ be a cancellable element of the lattice $\mathbb{SEM}$. If $\mathbf V$ satisfies some permutational identity of length $n\ge 3$ then it satisfies all identities of the form $\mathbf u\approx 0$ where $\mathbf u$ is a word of length $n$ depending on $<n$ letters.\qed
\end{corollary}

We say that a semigroup variety has a \emph{degree} $n$ if all nilsemigroups in this variety nilpotent of degree $\le n$ and $n$ is the least number with such a property.

\begin{proof}[Proof of Theorem~\textup{\ref{canc-permut-3}}. Necessity.]
 Let $\mathbf V$ be a semigroup variety that satisfies a permutational identity of length~3 and is a cancellable element of the lattice $\mathbb{SEM}$. Clearly, $\mathbf V$ is a modular element of $\mathbb{SEM}$. In view of Proposition~\ref{nec}, $\mathbf V=\mathbf{M\vee N}$ where $\mathbf M$ is one of the varieties $\mathbf T$ or $\mathbf{SL}$, while  $\mathbf N$ is a nil-variety. Proposition~\ref{alternative for cancellable} implies that $\mathbf V$ and therefore, $\mathbf N$ satisfies the identities~\eqref{xyz=yxz=xzy}. Now Corollary~\ref{non-linear = 0} applies with the conclusion that the identity~\eqref{xxy=0} holds in $\mathbf N$.

\medskip

\emph{Sufficiency.} Let $\mathbf V=\mathbf{M\vee N}$ where $\mathbf M$ is one of the  varieties $\mathbf T$ or $\mathbf{SL}$, while the variety $\mathbf N$ satisfies the identities~\eqref{xyz=yxz=xzy} and~\eqref{xxy=0}. We need to verify that $\mathbf V$ is cancellable. In view of Lemma~\ref{join with SL} we can assume that $\mathbf V=\mathbf N$. The identities~\eqref{xyz=yxz=xzy} imply $xyz\approx zyx$. Proposition~\ref{mod-permut-3} implies then that the variety $\mathbf N$ is a modular element of $\mathbb{SEM}$. Suppose that $\mathbf N$ is not a cancellable element of $\mathbb{SEM}$.  Let $\mathbf Y$ and $\mathbf Z$ be semigroup varieties such that $\mathbf N\vee\mathbf Y=\mathbf N\vee\mathbf Z$ and $\mathbf N\wedge\mathbf Y=\mathbf N\wedge\mathbf Z$. In view of Lemma~\ref{modular non-cancellable}, there is a variety $\mathbf N'$ such that
$$
\mathbf{N'\subseteq N},\,\mathbf{N'\vee Y}=\mathbf{N'\vee Z}=\mathbf{Y\vee Z}\text{ and }\mathbf{N'\wedge Y}=\mathbf{N'\wedge Z}.
$$
It is evident that the identities~\eqref{xyz=yxz=xzy} and~\eqref{xxy=0} implies
\begin{equation}
\label{xyx=yxx=0}
xyx\approx yx^2\approx 0.
\end{equation}
Being a subvariety of \textbf N, the variety $\mathbf N'$ satisfies the identities~\eqref{xyz=yxz=xzy},~\eqref{xxy=0} and~\eqref{xyx=yxx=0}. We aim to verify that $\mathbf Y=\mathbf Z$. This completes the proof because contradicts the claim that $\mathbf N$ is not a cancellable element of $\mathbb{SEM}$. By symmetry, it suffices to check that $\mathbf{Z\subseteq Y}$. Let $\mathbf{u\approx v}$ be an arbitrary identity that holds in $\mathbf Y$. It suffices to prove that this identity holds in $\mathbf Z$.

First of all we note that if the identity $\mathbf{u\approx v}$ holds in $\mathbf N'$ then it holds also in $\mathbf N'\vee\mathbf Y=\mathbf N'\vee\mathbf Z$  and therefore, in $\mathbf Z$. In particular, this allows us to assume that one of the words $\mathbf u$ and $\mathbf v$, say $\mathbf u$, does not equal to~0 in $\mathbf N'$. On the other hand, the fact that $\mathbf N'$ satisfies the identities~\eqref{xyz=yxz=xzy},~\eqref{xxy=0} and~\eqref{xyx=yxx=0} implies that every non-linear word except $x^2$ equals to~0 in $\mathbf N$ and therefore, in $\mathbf N'$. Thus, we can assume that either the word $\mathbf u$ is linear or $\mathbf u=x^2$. The latter case can be considered by literally repeating arguments from the consideration of the same case in the proof of sufficiency of Theorem~1.1 in~\cite{Gusev-Skokov-Vernikov-18+}. Suppose now that the word $\mathbf u$ is linear. The case when the identity $\mathbf{u\approx v}$ is not permutational also can be considered by literally repeating arguments from the consideration of the same case in the proof of sufficiency of Theorem~1.1 in~\cite{Gusev-Skokov-Vernikov-18+}. Finally, suppose that the identity $\mathbf{u\approx v}$ is permutational. If its length is greater than 2 then the identity $\mathbf{u\approx v}$ holds in $\mathbf N$ and therefore, in $\mathbf N'$. We reduce the proof to the case then the identity $\mathbf{u\approx v}$ is the commutative law.

Thus, the variety $\mathbf Y$ is commutative. The identity $xy\approx yx$ holds in the variety $\mathbf N'\wedge\mathbf Y=\mathbf N'\wedge\mathbf Z$. Therefore, there is a deduction of this identity from identities of the varieties $\mathbf N'$ and $\mathbf Z$. In particular, one of these varieties satisfies an identity of the form $xy\approx\mathbf w$. Suppose that this identity holds in $\mathbf N'$. Then Lemma~\ref{split} implies that $\mathbf N'$ either is commutative or satisfies the identity $xy\approx 0$. Clearly, $\mathbf N'$ is commutative always. Therefore, the commutative law holds in the variety $\mathbf N'\vee\mathbf Y=\mathbf N'\vee\mathbf Z$. Then the variety $\mathbf Z$ is commutative, and we are done.

It remains to consider the case when the identity $xy\approx\mathbf w$ holds in $\mathbf Z$. Lemma~\ref{split} implies that in this case either the variety $\mathbf Z$ is commutative or all nilsemigroups in $\mathbf Z$ are semigroups with zero multiplication. In the former case we are done. It remains to consider the latter case. In other words, $\mathbf Z$ is a variety of degree $\le 2$. Repeating literally arguments of the proof of~\cite[Lemma~3.1]{Gusev-Skokov-Vernikov-18+}, we can verify that the varieties $\mathbf Y$ and $\mathbf Z$ has the same degree. Therefore, $\mathbf Y$ is a variety of degree $\le 2$ too. By~\cite[Lemma~3]{Golubov-Sapir-82} or~\cite[Proposition~2.11]{Vernikov-08} this means that each of the varieties $\mathbf Y$ and $\mathbf Z$ saisfies one of the identities
\begin{align}
\label{left} xy&\approx x^{m+1}y,\\
\label{right} xy&\approx xy^{m+1},\\
\label{center} xy&\approx (xy)^{m+1}
\end{align}
for some natural $m$. Clearly, if $\mathbf Y$ [respectively $\mathbf Z$] satisfies one of these identities with $m=r$ [respectively $m=s$] then both the varieties satisfy corresponding identities with $m=rs$. Thus, we can assume that values of the parameter $m$ for $\mathbf Y$ and $\mathbf Z$ coincide. Further considerations are divided into seven cases.

\medskip

\emph{Case} 1: the variety $\mathbf Z$ satisfies the identity~\eqref{center}. Being commutative, the variety $\mathbf Y$ satisfies the identity $(xy)^{m+1}\approx(yx)^{m+1}$. Therefore, this identity holds in $\mathbf N'\vee\mathbf Y=\mathbf N'\vee\mathbf Z$. In particular, it holds in $\mathbf Z$. Hence $\mathbf Z$ satisfies the identities $xy\approx(xy)^{m+1}\approx(yx)^{m+1}\approx yx$.

\medskip

\emph{Case} 2: the varieties $\mathbf Y$ and $\mathbf Z$ satisfy the identity~\eqref{left}. Here $\mathbf Y$ satisfies the identities $x^{m+1}y\approx xy\approx yx\approx y^{m+1}x$. Therefore, the identity $x^{m+1}y\approx y^{m+1}x$ holds in $\mathbf N'\vee\mathbf Y=\mathbf N'\vee\mathbf Z$. In particular, it holds in $\mathbf Z$. Hence $\mathbf Z$ satisfies the identities $xy\approx x^{m+1}y\approx y^{m+1}x\approx yx$.

\medskip

\emph{Case} 3: the varieties $\mathbf Y$ and $\mathbf Z$ satisfy the identity~\eqref{right}. This case is dual to the previous one.

\medskip

\emph{Case} 4: the varieties $\mathbf Y$ and $\mathbf Z$ satisfy the identities~\eqref{left} and~\eqref{right} respectively. Here $\mathbf Y$ satisfies the identities $xy^{m+1}\approx y^{m+1}x\approx yx\approx xy\approx x^{m+1}y\approx yx^{m+1}$. Therefore, the identity $xy^{m+1}\approx yx^{m+1}$ holds in $\mathbf N'\vee\mathbf Y=\mathbf N'\vee\mathbf Z$. In particular, it  holds in $\mathbf Z$. Hence $\mathbf Z$ satisfies the identities $xy\approx xy^{m+1}\approx yx^{m+1}\approx yx$.

\medskip

\emph{Case} 5: the varieties $\mathbf Y$ and $\mathbf Z$ satisfy the identities~\eqref{right} and~\eqref{left} respectively. This case is dual to the previous one.

\medskip

\emph{Case} 6: the varieties $\mathbf Y$ and $\mathbf Z$ satisfy the identities~\eqref{center} and~\eqref{left} respectively. Here $\mathbf Y$ satisfies the identity $x^{m+1}y\approx(x^{m+1}y)^{m+1}$. Therefore, this identity holds in $\mathbf N'\vee\mathbf Y=\mathbf N'\vee\mathbf Z$ and therefore, in $\mathbf Z$. Hence $\mathbf Z$ satisfies the identities $xy\approx x^{m+1}y\approx(x^{m+1}y)^{m+1}\approx(xy)^{m+1}$. We reduce the situation to the Case 1 considered above.

\medskip

\emph{Case} 7: the varieties $\mathbf Y$ and $\mathbf Z$ satisfy the identities~\eqref{center} and~\eqref{right} respectively. This case is dual to the previous one.
\end{proof}

\medskip

config: arxiv


TITLE: Proof without words that $\tan ^{-1} 2-\tan ^{-1}\left(\frac{3}{2}\right)=\cot ^{-1} 8$
QUESTION [1 upvotes]: After reading the book “Proof without words” written by Mr Roger B. Nelsen, I started to collect questions with proof without words.  When I encounter problems now, I’ll first try to find the proof without words. For this one, I found one myself and want to share it with you now.

Is it beautiful?
Have you got any other alternate proofs?  Share with us if you have any.

REPLY [0 votes]: As suggested by Steven Stadnicki, I have just succeeded in completing a proof without words$$
\tan ^{-1} u-\tan ^{-1} v=\tan ^{-1}\left(\frac{u-v}{1+u v}\right)
$$
by the following diagram:


TITLE: A property of Radon-Nikodym derivatives
QUESTION [1 upvotes]: Is the Radon-Nikodym derivative linear in countably infinitely many input measures? That is, if $f_i$ is a Radon-Nikodym derivative of the measure $\nu_i$ w.r.t. $\mu$ for $i $ from $1$ to $ \infty $, then is it true that the sum of the $f_i$'s is a Radon-Nikodym derivative of the sum of the $\nu_i$'s w.r.t $\mu$?

REPLY [1 votes]: Yes. This is easy to see using the definitions and the monotone convergence theorem.
Simply note that
$$
\int_A \sum_i f_i d\mu = \sum_i \int_A f_i d\mu = \sum_i \nu_i(A) = (\sum_i \nu_i )(A)
$$
for each measurable $A$, where the first step is due to the monotone convergence theorem and the last step due to definition of $\sum_i \nu_i$.


text:
\begin{document}

\title{A Meshfree Generalized Finite Difference Method for Surface PDEs}

\runningheads{Suchde \textit{et~al.}}{A Meshfree GFDM for Surface PDEs}

\author{Pratik Suchde \affil{1}\corrauth, J\"org Kuhnert\affil{1}}

\address{\affilnum{1}Fraunhofer ITWM, 67663 Kaiserslautern, Germany}

\corraddr{P. Suchde. E-mail: pratik.suchde@itwm.fraunhofer.de}

\begin{abstract}
In this paper, we propose a novel meshfree Generalized Finite Difference Method~(GFDM) approach to discretize PDEs defined on manifolds. Derivative approximations for the same are done directly on the tangent space, in a manner that mimics the procedure followed in volume-based meshfree GFDMs. As a result, the proposed method not only does not require a mesh, it also does not require an explicit reconstruction of the manifold. In contrast to some existing methods, it avoids the complexities of dealing with a manifold metric, while also avoiding the need to solve a PDE in the embedding space. A major advantage of this method is that all developments in usual volume-based numerical methods can be directly ported over to surfaces using this framework. We propose discretizations of the surface gradient operator, the surface Laplacian and surface Diffusion operators. Possibilities to deal with anisotropic and discontinous surface properties (with large jumps) are also introduced, and a few practical applications are presented.
\end{abstract}

\keywords{Meshfree; Surface; Manifold; GFDM; Finite Difference}

\maketitle

\vspace{-6pt}

\section{Introduction}

The solution to partial differential equations~(PDEs) defined on a surface or manifold is of fundamental interest in various fields. They have application in the fields of computer graphics~\cite{Turk1991}, image processing~\cite{Diewald2000}, fluid flow~\cite{Myers2002}, and cell biology~\cite{Novak2007}, to name a few.

Most existing surface PDE solvers can be classified into two types. The first, referred to as intrinsic methods, solve the PDE directly on the manifold using either a mesh~\cite{Du2003,Du2006,Dziuk2013,Olshanskii2009}, a parametrization of the manifold~\cite{Floater2005}, or an explicit reconstruction of the manifold \cite{Liang2013}. In the second type of methods, referred to as embedding methods, the surface PDE is extended to a PDE defined on a band around the manifold, which is then discretized \cite{Bertalmio2001,Chu2018, Ratz2006,Ruuth2008}. The significant disadvantage of these methods is that they rely on discretizing a higher dimensional space, and thus they tend to be expensive in terms of computational time. On the other hand, the intrinsic methods have the disadvantage that the parametrization can be non-trivial to obtain, and that they need to deal with singularities arising in metric terms of surface differential operators~\cite{Fuselier2013}. A more detailed breakdown of existing literature on surface PDEs can be found in \cite{ChenThesis_CPM,IngridThesis_CPM}.

Constructing a good meshing of a manifold can be a very difficult process. It is thus often desirable to discretize a manifold directly with a set of scattered points, referred to as a point cloud. As a result, the need for meshfree methods to solve surface PDEs arises. In this paper, we present a new meshfree method for solving surface PDEs.

Several meshfree methods for surface PDEs have already been proposed. Most notable among these are the radial-basis function~(RBF) based methods \cite{Flyer2014,Fuselier2013}. RBF based surface PDE solvers combine the advantages of both intrinsic and embedding methods. They solve surface PDEs on the surface itself, without needing to parametrize or mesh the surface. However, they suffer from the drawback of needing a somewhat ad-hoc choice of basis function and related shape parameter, which affects stability and the conditioning of the linear systems \cite{Fuselier2013,Mongillo2011}. This choice is application dependent \cite{Simonenko2014}, given data~(right hand side or initial condition) dependent \cite{Davydov2011} , and can even be domain dependent. The optimal choice is not always known. We use a meshfree generalized finite difference method~(GFDM) to avoid this issue.

Meshfree GFDMs \cite{Fan2018,Gavete2017,Katz2010,Luo2016} are strong form meshfree methods that have been shown to be robust methods, and have been used in a wide variety of applications \cite{Drumm2008,Jefferies2015,Moller2007,Tramecon2013}. Approximations are carried out using a weighted least squares approach. In this paper, we propose a meshfree GFDM approach to solve PDEs defined on a surface. This is done by projecting local neighbourhoods to the tangent space, and performing approximations there. This extends the work of Demanet~\cite{Demanet2006} for mesh-based surface Laplacians to a meshfree context applicable to various differential operators. The proposed new method retains the advantages of RBF based methods, while avoiding the disadvantage of the ad-hoc choice of basis functions and shape parameters.

We note that meshfree GFDM approaches to solve surface PDEs have already been proposed by Liang \textit{et~al.} \cite{Liang2012,Liang2013}. In contrast, we do not rely on surface-based metrics, and thus avoid the issue of arising singularities. Another key difference from \cite{Liang2012,Liang2013} is the formulation used here enables transferring developments from volume-based numerical methods directly to surfaces.

The remainder of the paper is organized as follows. In Section~\ref{sec:Definitions}, we present the basics of meshfree GFDMs and introduce the notation used in the paper. Section~\ref{sec:SurfaceDOs} introduces a novel way to discrete surface differential operators in a meshfree GFDM setting. Section~\ref{sec:BoundaryConditions} contains a short note on the implementation of boundary conditions. While Section~\ref{sec:SurfaceDOs} deals with $2$-manifolds in $3$ dimensional space,  Section~\ref{sec:HighD_coD} extends those ideas to higher dimensions and co-dimensions. Section~\ref{sec:NumResults} presents a range of numerical examples and validation, and the paper is concluded with a discussion on the work in Section~\ref{sec:Conclusion}.

\section{Preliminaries}
\label{sec:Definitions}

\subsection{Basics and Nomenclature}

For the majority of this paper, we consider a smooth orientable $2$-manifold $M$ embedded in $\mathbb{R}^3$. Extensions of the ideas presented here to higher dimensions and co-dimensions can be done easily, and will be discussed briefly in Section \ref{sec:HighD_coD}. The manifold $M$ is also referred to as the embedded space or surface, while $\mathbb{R}^3$ is also referred to as the embedding space.

We consider a manifold given by an unevenly distributed point cloud consisting of $N$ points.  For each point $i$ on the manifold, approximations are carried out on the neighbourhood or support $S_i$ consisting of $|S_i|$ nearby points. $S_i$ is based on proximity, given by Eucledian distances in the embedding space. Throughout this paper, all distances are computed in the embedding space only. Distances along the manifold are never needed, and are not computed. The size of the support $S_i$ is given by the smoothing length $h_i$. We adopt the following distance conventions from volumetric meshfree GFDM for fluid flow \cite{Drumm2008, Suchde2017_CCC}. During set up of the point cloud, it is ensured that no two points are closer than $r_{min}h$, and that there is no hole of size $r_{max}h$ on the surface that does not contain any point. The parameters $r_{min}=0.2$ and $r_{max}=0.45$ are fixed across all simulations, and are adopted from volumetric meshfree GFDM conventions \cite{Suchde2017_CCC}. This ensures that that each support of size $h$ has about $15-20$ points. This also results in $h$ serving as an indication of the point cloud spacing. In all simulations, it is assumed that $h$ is chosen such that the local neighbourhoods are sufficiently dense to accurately capture the local curvature of the surface.

We assume that at each interior point $i$, the unit surface normal $\vec{n}_i$ and unit tangents $\vec{t}_{1,i}$, $\vec{t}_{2,i}$ to the surface have already been computed. Similarly, for each boundary point, the unit surface normal $\vec{n}_i$, the unit surface tangent $\vec{t}_i$, and the unit boundary normal $\vec{\nu}_i$ are assumed to be known. Note that these normals and tangents form an orthogonal system of vectors. Normal and tangent computation could have been done in any of multiple ways. One possibility is to use Principle Component Analysis~(PCA) to construct the normal and tangent information based on the eigenvalues of local covariance matrices (see, for example, \cite{Liang2013,Mitra2003}). Weighted PCA approaches \cite{Petronetto2013} could also be used for the same. Alternatively, surface normals and tangents could be computed based on local geometric information in each neighbhourhood using the same procedures done in meshfree GFDM simulations of fluids to compute normals at the free surface (see, for example, \cite{Edgar2017}). This latter approach will be used for all simulations in this paper. Further, if the point cloud is based directly on a CAD model, normal information could also be taken from the CAD model itself, if it is available, and the tangents could be computed accordingly.

The following notation is used for differential operators throughout this paper. $\nabla$, $\Delta$ are used to denote the continuous operators in the embedding space $\mathbb{R}^3$. Subscripts $M$ and $T$ to the differential operators are used to indicate the corresponding operator defined on the manifold and the tangent space respectively, while the subscript $tn$ is used to indicate operators in $\mathbb{R}^3$ in a rotated tangent-normal frame. A $\tilde{}\,$ is used overhead to explicitly state, if needed, that the operator being considered is a discrete one. A second subscript $i$ is used to indicate that the discrete operator is being considered at the point $i$. Thus, $\widetilde{\Delta}_{M,i}$ indicates the discrete surface Laplacian at point $i$.

\subsection{Volume-based Meshfree GFDMs}
\label{sec:VolumetricGFDMs}

Classical meshfree GFDMs for volume or bulk domains shall henceforth be referred to as `volumetric' GFDMs. They are strong-form methods, in which for a function $u$ defined at each numerical point $i=1,2,\dots,N$, its derivatives are approximated as

\begin{align}
	\nabla u(\vec x_i) &\approx \widetilde{\nabla}_i u =
	\left( \def\arraystretch{2.0} \begin{array}{c}
		\sum_{j\in S_i}c_{ij}^x u_j	\\
		\sum_{j\in S_i}c_{ij}^y u_j	\\
		\sum_{j\in S_i}c_{ij}^z u_j	\\
	\end{array}  \right) \,,\label{Eq:GFDM_Definition1}\\

	\Delta u(\vec x_i) &\approx \widetilde{\Delta}_i u = \sum_{j\in S_i}c_{ij}^\Delta u_j\,,\label{Eq:GFDM_Definition2}
\end{align}

where $\widetilde{(\cdot)}_i$ indicates the discrete differential operator at point $i$, and the differential operators without the tilde indicate the continous differential operators. All meshfree GFDM notation used here follows from \cite{Suchde2017_CCC}. For each point $i$, the stencil coefficients $c_{ij}$ are found using a weighted least squares approach. The weighted sum of the stencil coefficients is minimized such that monomials $m\in\mathcal{P}$ up to a certain order, usually $2$, are exactly differentiated. For example, for the Laplacian,

\begin{align}
	\sum_{j\in S_i}c_{ij}^{\Delta}m_j &= \Delta m (\vec{x_i}) \qquad \forall m\in\mathcal{P}\,,\label{Eq:Consistency}\\
	\text{min } J_i &= \sum_{j\in S_i} \left( \frac{c_{ij}^{\Delta}}{W_{ij}} \right)^2\,, \label{Eq:BasicMin}
\end{align}

where $W$ is a weighting function. Throughout this paper, we use a Gaussian weighting function

\begin{equation}
	W_{ij}=\exp(-W_F \frac{\|\vec{x}_j-\vec{x}_i\|^2}{h_i^2 + h_j^2})\,,
\end{equation}

where $W_F>0$, and the inclusion of the smoothing length $h_j$ of the neighbour point is for point clouds with uneven point distributions. Note that for a central point $i$, the weighting function is only relevant for the neighbouring points $j\in S_i$. The same procedure also holds for other differential operators. The monomials depend on distance relative to the central point $i$. So, for example, in $\mathbb{R}^2$, for the $x$ derivative, the monomial consistency conditions lead to the following system

\begin{equation}
\label{Eq:DO_Example}
\underbrace{\left(\begin{array}{ccc}
	\cdots& 1 &\cdots\\
	\cdots& \delta x_{ij} &\cdots\\
	\cdots& \delta y_{ij} &\cdots\\
	\cdots& \delta x_{ij}^2 &\cdots\\
	\cdots& \delta y_{ij}^2 &\cdots\\
	\cdots& \delta x_{ij}\delta y_{ij} &\cdots\\
	\end{array}\right)}_{M_i^T}

\underbrace{\left(\begin{array}{c}
	\vdots \\
	c_{ij}^x\\
	\vdots\\
	\end{array}\right)}_{\vec{c}_i^{\,x}} =

\underbrace{\left(\begin{array}{c}
	0\\
	1\\
	0\\
	0\\
	0\\
	0\\
	\end{array}\right)}_{\vec{b}_i^{\,x}}\,,
\end{equation}

where $\delta x_{ij}=x_j - x_i$ and $\delta y_{ij} = y_j - y_i$. The number of neighbours is always taken to be larger than the number of monomial consitency conditions. Continuing the example of the $x$ derivative, the minimization equivalent of Eq.\,\eqref{Eq:BasicMin} can be written as $\text{min } J_i = \|W_i^{-1}\vec{c}_i^{\,x}\|^2$, with $W_i$ being the diagonal matrix of weights. This minimization leads to the differential operators given by $\vec{c}_i^{\,x} = W_i^2M_i\left(M_i^TW_i^2M_i\right)^{-1}\vec{b}_i^{\,x}$. Note that the LHS matrix in Eq.\,\eqref{Eq:DO_Example} is the same for all differential operators being approximated. Thus, multiple differential operators can be computed numerically with a single minimization procedure with multiple right hand sides~\cite[Section A.2]{Suchde2018_Thesis}.

We wish to use a similar approach here for surface derivatives. All numerical differential operators on a manifold will be defined in a manner similar to Eq.\,\eqref{Eq:GFDM_Definition1} and Eq.\,\eqref{Eq:GFDM_Definition2}, and the procedure for computing the stencil coefficients will be done in a similar manner to Eq.\,\eqref{Eq:Consistency} and Eq.\,\eqref{Eq:BasicMin}.

\section{Surface Differential Operators using a Meshfree GFDM}
\label{sec:SurfaceDOs}

Consider a function $u: M \rightarrow \mathbb{R}$ defined on the surface, and an open subset $\Omega$ of $\mathbb{R}^3$ containing $M$ i.e. $M \subset \Omega \subset \mathbb{R}^3$. A function $\hat{u}: \Omega \rightarrow \mathbb{R}$ is said to be an extension of $u$ if $u$ and $\hat{u}$ agree on the manifold, $\hat{u}|_{M} = u$. Further, $\hat{u}$ is said to be a normal extension of $u$ if $\vec{n}\cdot\nabla \hat{u} \equiv 0$.

In this paper, we use normal extensions to extend functions \textit{locally}. We then use these extensions to define numerical differential operators entirely on the tangential plane at each point. We emphasize that the extension is \textit{not} done numerically as is the case for embedding methods such as the closest point method~\cite{Marz2012,Ruuth2008}. Further, no dummy or virtual points are needed for differential operator compuation, as is the case for some RBF-based methods for surface PDEs~\cite{Piret2012}. Here, only the manifold $M$ is discretized, whereas $\Omega$ is not. The concept of (normal) extensions is only used to derive a way to compute differential operators directly on the tangential plane in a straight forward manner.

\subsection{Surface Gradient}
\label{sec:SurfaceGradient}

The surface gradient of a function can be defined as the conventional gradient of its extension with the component normal to the surface removed

\begin{align}
	\nabla_{M} u &= \nabla \hat{u} - \vec{n} \left( \vec{n} \cdot \nabla \hat{u} \right)\,, \label{Eq:SurfaceGradient_Basic}\\
	&= \left( \mathbf{P}\nabla \right) \hat{u}\,,
\end{align}

where $\mathbf{P}=\mathbf{I} - \vec{n}\vec{n}^T$ for identity $\mathbf{I}$ is the projection matrix (to the tangent space). Here, we define $\hat{u}$ such that it is a normal extension of $u$. Thus, $\vec{n}\cdot\nabla \hat{u}= 0$, and Eq.\,\eqref{Eq:SurfaceGradient_Basic} reduces to

\begin{equation}
	\label{Eq:GradientTheorem}
	\nabla_{M} u = \nabla \hat{u} \,.
\end{equation}

Thus, determining a numerical approximation $\widetilde{\nabla}\hat{u}$ to $\nabla\hat{u}$  (the gradient of the extension $\hat{u}$ ),  will give us an approximation to the surface gradient $\widetilde{\nabla}_{M}u$ of the original function $u$.

To compute the approximation $\widetilde{\nabla}\hat{u}$, we once again make use of the fact that $\vec{n}\cdot\nabla \hat{u}= 0$ everywhere on the manifold. Rather than computing $\widetilde{\nabla}\hat{u}$ directly, we first consider the rotated components

\begin{equation}
	\nabla_{tn}\hat{u} = \left(\begin{array}{c}
	\vec{t}_1\cdot\nabla \hat{u} \\
	\vec{t}_2\cdot\nabla \hat{u} \\
	\vec{n}\cdot\nabla \hat{u} \\
	\end{array}\right) = \left(\begin{array}{c}
	\vec{t}_1\cdot\nabla \hat{u} \\
	\vec{t}_2\cdot\nabla \hat{u} \\
	0 \\
	\end{array}\right) \,,
\end{equation}

where the subscript $tn$ indicates the rotated $\vec{t}_1$, $\vec{t}_2$, $\vec{n}$ coordinate frame at a particular point. Thus, we only need to compute numerical approximations to the tangential components $\vec{t}_k \cdot \nabla \hat{u}$ for $k=1,2$. Once these tangential components are known, they can be rotated to get the gradient $\nabla \hat{u}$.

\begin{equation}
	\label{Eq:RotationOfGradient}
	\nabla \hat{u} =  \underbrace{ \left(\begin{array}{ccc}
	\vec{t}_1 & \vec{t}_2 & \vec{n}
	\end{array}\right)}_{R^T}

	\nabla_{tn}\hat{u} \,,
\end{equation}

where the rotation matrix $R^T$ is composed of the the tangential and normal column vectors.

\subsubsection{Numerical Surface Gradient Operator}
~\\
We showed above that the problem of computing numerical differential operators for the surface gradient of a function can be reduced to that of computing numerical differential operators for the tangential components of the regular volumetric gradient of the normal extension of the function.

Numerically, at each discrete point $i$ on the manifold, we compute approximations of the tangential components of $\nabla \hat{u}$ entirely on the tangent plane $T_i$ spanned by $\vec{t}_{1,i}$ and $\vec{t}_{2,i}$, i.e. we compute

\begin{equation}
	\widetilde{\nabla}_{T} u = \left(\begin{array}{c}
	\vec{t}_1\cdot \widetilde{\nabla} \hat{u} \\
	\vec{t}_2\cdot \widetilde{\nabla} \hat{u}
	\end{array}\right) \,,
\end{equation}

where $\widetilde{\nabla}_{T}$ represents the discrete $2$ dimensional gradient in the tangential plane. For this, we project each point $j \in S_i$ to the tangential plane $T_i$. Let the projection of the point $j \in S_i$ to $T_i$ be labelled as $j_{T_i}$. Since $u$ is being extended normally (the derivative in the normal direction is $0$), $\hat{u}$ evaluated at point $j_{T_i}$  is the same as $u$ evaluated at point $j$. Thus, we set

\begin{equation}
	\label{Eq:Projection_j_Equiv}
	\hat{u}_{j_{T_i}} = u_j
\end{equation}

A neighbouring point $j \in S_i$ is projected to the tangent plane $T_i$ along the surface normal of the central point $\vec{n}_i$. As a result, tangential distances are preserved. The distances between the central point $i$ and the projected points $j_{T_i}$ can be computed by simply rotating the original distances. If $\delta\vec{x}_{ij} = \vec{x}_j - \vec{x}_i$, then the distances in the tangential and normal directions are given simply by $R\, \delta\vec{x}_{ij}$, where $R$ is the rotation matrix introduced in Eq.\,\eqref{Eq:RotationOfGradient}. Of these, only the tangential distances are needed for derivative computation. The computation of the differential operators according to a procedure like Eq.\,\eqref{Eq:Consistency} and Eq.\,\eqref{Eq:BasicMin} only requires the distances between points, as shown in Eq.\,\eqref{Eq:DO_Example}. Since we already have these distances, there is no need to actually add a point numerically at the projected location.

We note that the idea of projecting to the tangent space to compute numerical approximations is not a novel one. It has also been done by mesh-based surface PDE solvers. Lai \textit{et~al.} \cite{Lai2013} do the same, but while also working with surface-based metrics. Demanet~\cite{Demanet2006} does the same in the mesh-based framework, but only for the surface Laplacian. There, they project points along the surface normal of the neighbouring point itself $\vec{n}_j$. Such a procedure would involve a little more work numerically, and as we show later, is less accurate than the central normal projection used here. The difference between the two methods of projecting neighbouring points to the tangent space is illustrated in Figure~\ref{Fig:Projection} for a $1$-dimensional manifold in $\mathbb{R}^2$.

Unlike the neighbour normal projection, the central normal projection method does not give a global normal extension of the function. i.e. if the union of all the projected locations on tangential planes $\mathcal{T} = \cup_{i,j\in S_i}\, \vec{x}_{j_{T_i}}$ is considered,
$\hat{u}$ defined on $\mathcal{T}$ does not form a global normal extension of $u$. However, this is not relevant in the present context. Our interest is only local. For each point $i$, a virtual projection of its neighbours in $S_i$ to its tangential plane $T_i$ is only performed to numerically compute the derivatives at point $i$ itself.

\begin{figure}
  \centering

  \includegraphics[width=0.45\textwidth]{./Projection1}
  \includegraphics[width=0.45\textwidth]{./Projection2}

  \caption{Projection of neighbouring points to the tangent space along the central normal~(left), and the neighbour normals~(right). The central point is shown with an additional circle around it, and all its neighbouring points are marked in red. The remaining black points on the manifold are shown for reference. The projected locations on the tangent space are marked in blue. The manifold is shown in black, while the tangent line at the central point is shown with a blue dashed line.}

  \label{Fig:Projection}
\end{figure}

For each point $i$, once the tangential distances to the neighbouring points are known, we  compute $2$ dimensional volumetric numerical differential operators on the tangential plane for the first derivatives along the $\vec{t}_{1,i}$ and $\vec{t}_{2,i}$ directions, as explained in Section~\ref{sec:VolumetricGFDMs}.

\begin{align}
	\sum_{j\in S_i}c_{ij_{T}}^{t_k}m_{j_{T}} &= \frac{\partial}{\partial t_k} m (\vec{x}_i) \qquad \forall m\in\mathcal{P}_{T}\,,\label{Eq:TP_Consistency}\\
	\text{min } J_i &= \sum_{j\in S_i} \left( \frac{c_{ij_{T}}^{t_k}}{W_{ij_{T}}} \right)^2\,, \label{Eq:TP_Min}
\end{align}

for $k=1,2$, where $\mathcal{P}_{T}$ are the set of monomials, usually up to order $2$, in $\vec{t}_{1,i}$ and $\vec{t}_{2,i}$ on the tangent plane.

Using Eq.\,\eqref{Eq:RotationOfGradient} and Eq.\,\eqref{Eq:GradientTheorem}, the computed differential operator stencil coefficients can now be rotated to obtain the numerical surface gradient operator. We have

\begin{equation}
	\label{Eq:SurfGrad_FullDefinition}
	\widetilde\nabla_{M,i} u =
									\left(\def\arraystretch{2.0} \begin{array}{c}
									\sum_{j \in S_i}c_{ij}^{M,x} u_j \\												\sum_{j \in S_i}c_{ij}^{M,y} u_j \\
									\sum_{j \in S_i}c_{ij}^{M,z} u_j \\
										\end{array}\right)\,,
\end{equation}

where $c_{ij}^{M,x}$ are the stencil coefficients for the surface gradient in the $x$ direction, and similarly for the other directions. Further,

\begin{equation}
	\left(\def\arraystretch{2.0} \begin{array}{c}
	c_{ij}^{M,x} \\														   			c_{ij}^{M,y} \\
	c_{ij}^{M,z} \\
		\end{array}\right) =

	R^T 		\left(\def\arraystretch{2.0} \begin{array}{c}
	c_{ij}^{t_1} \\														   						c_{ij}^{t_2} \\
	c_{ij}^{n} \\
		\end{array}\right)\,,
\end{equation}

for $c_{ij}^n = 0$.

This procedure used for the computation of the numerical surface gradient operators can be extended easily to obtain discretizations for any surface differential operator.

A key point to note here is that the main computation of differential operator stencil coefficients, according to Eq.\,\eqref{Eq:TP_Consistency} and Eq.\,\eqref{Eq:TP_Min}, follows the exact same procedure as that of regular volumetric differential operator computation. As a result, any and all developements in volumetric meshfree GFDMs can be directly used for surface PDEs. This is the one of the biggest advantages of this method over existing work for the same. A significant amount of work has been done to modify GFDM differential operator for different ends, and they can be easily carried over to surface PDEs using the present work. These include developments such as higher order spatial discretizations~\cite{Milewski2012}, conservation~\cite{Chiu2012, Suchde2017_CCC}, accuracy considerations~\cite{Suchde2018_INSE}, upwinding methods for advection \cite{Praveen2007, SeifarthThesis}, staggered methods \cite{Trask2017}, among others. A few examples of carrying over volumetric GFDMs developments to surface PDEs are shown in the coming sections.

\subsection{Surface Divergence}
\label{sec:SurfaceDivergence}

Consider a vector valued function $\vec{v} = (v^1, v^2, v^3)$ defined on the manifold. A normal extension of $\vec{v}$ is obtained by a normal extension of each component of $\vec{v}$. i.e. $\vec{\hat{v}} = \left( \hat{v}^1, \hat{v}^2, \hat{v}^3 \right)$. Now, the surface divergence of $\vec{v}$ can be written as

\begin{align}
	\nabla_{M}\cdot \vec{v} &=  \left(\mathbf{P}\nabla\right) \cdot \vec{\hat{v}}\,,\\
	&= \nabla\cdot\vec{\hat{v}} - \vec{n} \cdot \left[ \left( \vec{n} \cdot \nabla \right) \vec{\hat{v}} \right] \,,\\
	&= \nabla\cdot\vec{\hat{v}} \label{Eq:DivergenceTheorem}\,.
\end{align}

The numerical gradient operators computed above can be used to compute the surface divergence. For brevity, we write the numerical surface gradient according to Eq.\,\eqref{Eq:SurfGrad_FullDefinition} as

\begin{equation}
	\label{Eq:SurfaceGradientConcise}
	\widetilde{\nabla}_{M,i}u = ( G_{1i} u, G_{2i} u, G_{3i} u)^T\,,
\end{equation}

where $G_{ki}$ are the operators computed in the previous section. For example, $G_{1i}u = \sum_{j\in S_i}c_{ij}^{M,x}u_j$. Using the notation of Eq.\,\eqref{Eq:SurfaceGradientConcise}, the numerical surface divergence can be computed as

\begin{equation}
	\widetilde{\nabla}_{M}\cdot \vec{v} = \sum_{k=1,2,3} G_{ki} v^k\,.
\end{equation}

\subsection{Surface Laplacian}
\label{sec:SurfaceLaplace}

The surface Laplacian, or Laplace Beltrami, of a scalar valued function is defined as

\begin{equation}
	\label{Eq:LapBeltDefinition}
	\Delta_{M} u = \nabla_{M}\cdot \nabla_{M} u\,.
\end{equation}

Using Eq.\,\eqref{Eq:GradientTheorem} and Eq.\,\eqref{Eq:DivergenceTheorem} leads to

\begin{equation}
	\label{Eq:LaplaceTheorem}
	\Delta_{M} u = \Delta \hat{u}	\,.
\end{equation}

Similar to the surface gradient case, we compute a numerical approximation to $\Delta \hat{u}$ on the tangent plane. Using the same ideas as that in Section~\ref{sec:SurfaceGradient}, and specifically Eq.\,\eqref{Eq:RotationOfGradient}, we get

\begin{align}
	\Delta_{M} u &= \Delta \hat{u}\,,	\\
	&= \nabla \cdot \nabla \hat{u} \,, \\
	&= \nabla \cdot \left( R^T \nabla_{tn} \hat{u} \right)\,,\\
	&= \left( R \nabla \right) \cdot \left( R R^T \nabla_{tn} \hat{u} \right)\,, \label{Eq:DivRotInv1}\\
	&= \nabla_{tn} \cdot \left( R R^T \nabla_{tn} \hat{u} \right)\,,\\
	&= \Delta_{tn} \hat{u}\,,\\
	&= \Delta_{T} \hat{u}\,,
\end{align}

where $\Delta_{T}$ is the $2$ dimensional Laplacian on the tangent plane, and $RR^T = R^TR = I$. Eq.\,\eqref{Eq:DivRotInv1} uses the fact that the divergence is invariant under rotations. The last equation arises due to the $\vec{n}$ direction derivatives being $0$. The same can alternatively be derived using the rotational invariance of the Laplacian.

Similar to before, we compute a $2$ dimensional volumetric numerical operator for the Laplacian on the tangential plane, which directly gives the discretization for the surface Laplacian.

\begin{align}
	\sum_{j\in S_i}c_{ij_{T}}^{\Delta_{T}}m_{j_{T}} &= \Delta_{T} m (\vec{x}_i) \qquad \forall m\in\mathcal{P}_{T}\,,\label{Eq:TPL_Consistency}\\
	\text{min } J_i &= \sum_{j\in S_i} \left( \frac{ c_{ij_{T}}^{\Delta_{T}} } {W_{ij_{T}} } \right)^2\,. \label{Eq:TPL_Min}
\end{align}

The numerical surface Laplacian is then given by

\begin{equation}
	\label{Eq:LapBeltNumericalDefinition}
	\widetilde{\Delta}_{M} u = \sum_{j\in S_i} c_{ij}^{\Delta_M}u_j\,,
\end{equation}

where $c_{ij}^{\Delta_M} = c_{ij}^{\Delta_T}\;$.

\subsubsection{Optimized Surface Laplacian}
~\\
As mentioned earlier, all modifications of volume based meshfree GFDMs can be directly carried over to the case of surface operators. An important issue in volumetric meshfree GFDMs is the use of optimized Laplacian stencils to improve stability for Poisson problems \cite{SeiboldThesis}. For this, instead of a direct minimization according to Eq.\,\eqref{Eq:TPL_Min}, it is desired that the central stencil value $| c_{ii}^{\Delta} |$ is made ``as large as possible'' in relation to the neighbouring ones $| c_{ij}^{\Delta} |, \phantom{s} j\neq i$, while maintaining the consistency conditions of Eq.\,\eqref{Eq:TPL_Consistency}. Several procedures for the same have been done in the past (for example, \cite{Froese2018,Seibold2008}). Here, we follow the method done in our earlier work \cite[Section~2.5.5]{Suchde2018_Thesis}. A short explanation of the same is given in Appendix~\ref{App:Laplace}.

\subsection{Anisotropic Surface Laplacian}

We generalize the previous section to an anisotropic surface Laplacian~(surface diffusion) operator $\nabla_{M} \cdot \kappa \nabla_{M}\;$, for diffusion coefficeint $\kappa$. Proceeding in the same manner as earlier, the surface diffusion of a scalar valued function can be evaluated as follows

\begin{align}
	\nabla_{M} \cdot \left( \kappa \nabla_{M} u \right) &= \nabla\cdot\left( \hat{\kappa} \nabla \hat{u} \right)\,,\\
    & = \nabla \cdot \left(\hat{\kappa} R^T \nabla_{tn} \hat{u} \right)	\,,\\
    & = \left( R\nabla \right) \cdot \left( R\hat{\kappa} R^T \nabla_{tn} \hat{u} \right)	\,,\label{Eq:DivRotInv2}\\
    &= \nabla_{tn} \cdot \left( \underbrace{R\hat{\kappa} R^T}_{\hat{\kappa}_R  } \nabla_{tn} \hat{u} \right)\,,\\
	&= \nabla_{T} \cdot \left( \hat{\kappa}_{RT}  \nabla_{T} \hat{u} \right) \,,
\end{align}

where $\nabla_{T} \cdot \hat{\kappa}_{RT}  \nabla_{T}$ is the $2$ dimensional diffusion operator on the tangential plane, $\hat{\kappa}$ is the (possibly component-wise) normal extension of $\kappa$, and $\hat{\kappa}_{RT}$ is the appropriate submatrix of $\hat{\kappa}_{R}$. For scalar valued $\kappa$, $\hat{\kappa}_R = \hat{\kappa} = \hat{\kappa}_{RT}$. Eq.\,\eqref{Eq:DivRotInv2} uses the fact that the divergence is invariant under rotations. The last equation arises due to the $\vec{n}$ direction derivatives being $0$.

Thus, similar to the case in the earlier sections, the surface diffusion operator at a point $i$ has been reduced to the $2$ dimensional volumetric diffusion operator on the tangent plane of that point, which is approximated numerically.

\begin{align}
	\sum_{j\in S_i}c_{ij_{T}}^{D_{T}}m_{j_{T}} &= D_{T} m (\vec{x}_i) \qquad \forall m\in\mathcal{P}_{T}\,,\label{Eq:TPD_Consistency}\\
	\text{min } J_i &= \sum_{j\in S_i} \left( \frac{ c_{ij_{T}}^{D_{T}} } {W_{ij_{T}} } \right)^2\,, \label{Eq:TPD_Min}
\end{align}

where $D_{T}$ is used as a shorthand for $\nabla_{T} \cdot  \hat{\kappa}_{RT}  \nabla_{T}\;$ . The numerical surface diffusion operator is then given by

\begin{equation}
	\label{Eq:DiffusionNumericalDefinition}
	\nabla_{M} \cdot \left( \kappa \nabla_{M} u \right) \approx
	\widetilde{D}_{M} u = \sum_{j\in S_i} c_{ij}^{D_M}u_j\,,
\end{equation}

where $c_{ij}^{D_M} = c_{ij}^{D_T} \;$.

\subsubsection{Anisotropic Surface Laplacian with Large Jumps in Diffusion Coefficient}
\label{sec:DiffusionLargeJumps}
~\\
It is often required to model diffusion with not just a discontinuous diffusion coefficient, but one with large jumps~(with several orders of magnitude). To do the same, modifications needs to be made to the standard procedure of computing the numerical diffusion operator mentioned above.

To achieve the same on a manifold, we extend our earlier volumetric work. For a scalar-valued $\kappa$, the following consistency conditions are enforced in addition to Eq.\,\eqref{Eq:TPD_Consistency}.

\begin{align}
	\sum_{j\in S_i}c_{ij_{T}}^{D_{T}} \frac{1 }{\kappa_{j_{T} }} &= -\Delta_{M} \left( \log \kappa \right) \,,\label{Eq:DiffExtra1}\\
	\sum_{j\in S_i}c_{ij_{T}}^{D_{T}} \frac{\delta s_{ij_{T}} }{\kappa_{j_{T} }} &= - \frac{\partial}{\partial s_{M} } \left( \log\kappa \right) \,,\label{Eq:DiffExtra2}\\
		\sum_{j\in S_i}c_{ij_{T}}^{D_{T}} \frac{\left( \delta  s_{ij_{T}} \right)^2}{\kappa_{j_{T} }} &= 2 \,,\label{Eq:DiffExtra3}
\end{align}

where $\delta s_{ij_{T}}$ are distances along the direction of greatest change in $\kappa$, $\vec{s} = \frac{ \nabla_{M} \kappa}{ \| \nabla_{M} \kappa \|}$; and $\frac{\partial}{\partial s_{M} }$ denotes the directional surface derivative along $\vec{s}$, which is obtained numerically by rotating $\nabla_{M}\;$. The proof of the validity of the these extra conditions is given in Appendix~\ref{App:Diffusion}.

We note that this is similar in concept to the work of Yoon and Song \cite{Yoon2014} who add step functions, wedge functions and scissor functions to the polynomial test functions for GFDMs.

\section{Boundary Conditions}
\label{sec:BoundaryConditions}

One of the biggest advantages of meshfree GFDMs over particle-based meshfree methods such as SPH is the ease of handling a vast variety of boundary conditions \cite{Suchde2018_Thesis}. This is carried over to meshfree GFDMs on manifolds as well. To illustrate the same, we consider a surface Poisson equation

\begin{align}
	\Delta_M u &= f \qquad \text{in } M\,,\\
	u & = g \qquad \text{on } \partial M_1\,,\\
	\vec{\nu}\cdot\nabla_M u &= l \qquad \text{on } \partial M_2\,,
\end{align}

where $\partial M_1$ and $\partial M_2$ are parts of the manifold boundary. The discretized linear system would simply be given by

\begin{align}
	\sum_{j \in S_i} c_{ij}^{\Delta_M} u_j &= f_i \qquad \text{for } i \in M \setminus \partial M  \,,\\
	u_i &= g_i  \qquad \text{for } i \in \partial M_1 \,,\\
	\sum_{j \in S_i} c_{ij}^{M,\nu} u_j &= l_i \qquad \text{for } i \in  \partial M_2\,,
\end{align}

where $c_{ij}^{M, \nu} = \vec{\nu} \cdot \left( c_{ij}^{M,x}, c_{ij}^{M,y}, c_{ij}^{M,z}     \right)^T$.

We note that the neighbourhoods $S_i$ for boundary points $i$ will always be ``one-sided". Similar to the case of volumetric meshfree GFDMs, the computation of differential operators on boundary points follows the exact same procedure as that for interior points. The use of ghost or virtual nodes outside the domain is \textit{not} done. More details about implementing different kinds of boundary condition using meshfree GFDMs can be found in our earlier work \cite{Suchde2018_Thesis}.

\section{Higher Dimensions and Co-Dimensions}
\label{sec:HighD_coD}

The ideas presented in this paper can easily be extended to manifolds in higher dimensions, or higher co-dimensions. The main difference is the change of the projection operator. For a $k$ dimensional manifold in $\mathbb{R}^n$, we have

\begin{equation}
 	\mathbf{P} = \mathbf{I} - \sum_{r=1}^{n-k}\vec{n}_r\vec{n}_r^T\,,
\end{equation}

where the normal space is spanned by the unit normals $\vec{n}_r$. At the discrete level, once again, the only change is the way of projecting neighbouring points to the tangent space. Distances in the tangential space can still be computed by rotating actual distances in the embedding space. The only difference would be in the rotation matrix, which would computed as

\begin{equation}
	R^T = \left(\begin{array}{cccccc}
	\vec{t}_1 & \cdots & \vec{t}_k & \vec{n}_1 & \cdots & \vec{n}_{n-k}
	\end{array}\right)\,,
\end{equation}

where $\vec{t}_1, \dots ,\vec{t}_k$ are orthogonal to each other and each $\vec{n}_r$. Thus, $R^T$ has mutually orthogonal columns. Once the distances in the tangential space are known~(given by the appropriate submatrix of $R\, \delta\vec{x}_{ij} $), the differential operators can be computed in the same manner as volumetric $k$-dimensional operators.

\section{Numerical Results and Validation}
\label{sec:NumResults}

We emphasize that we do not use uniformly spaced point clouds, with the exception of the first example. All irregularly spaced point clouds are setup in a manner similar to that done in several volumetric meshfree GFDMs. Starting from a CAD file for the geometry, points are placed using an advancing front technique for point clouds, like that done by Drumm \textit{et~al.} \cite{Drumm2008}. Using the distance conventions used in \cite{Drumm2008, Suchde2017_CCC} results in about $15-20$ points in each neighbourhood. Normal computation follows the procedure used for free surfaces of volumetric flow by meshfree GFDM, as done by \cite{Edgar2017}.

Unless specified otherwise, in all numerical examples monomials up to the second order are considered in the computation of all numerical differential operators. Further, the sparse linear systems arising in each example are solved with a BiCGSTAB iterative solver \cite{BiCGSTAB} without the use of any pre-conditioner. If available, the solution of the previous time level is used as an initial guess.

\subsection{Heat Equation on a Sphere}

As a validation case, we consider the surface diffusion equation

\begin{equation}
	\label{Eq:HeatEquation}
	u_t = \Delta_{M} u \,,
\end{equation}

on a unit sphere. For initial conditions $u_0 = xy$, the analytical solution to Eq.\,\eqref{Eq:HeatEquation} is given by

\begin{equation}
	\label{Eq:AnalSurfHeat}
	u_{exact}(\vec{x},t) = \exp(-6t)xy \,.
\end{equation}

We note that the same example has also been considered in Chen \textit{et~al.} \cite{Chen2014}. The sphere is discretized with a quasi-uniformly distributed point cloud generated by DistMesh \cite{Persson2004}. Irregularly spaced point clouds are considered in the later sections. For consistency with the case of irregularly spaced point clouds considered from the next example, here, neighbourhoods are taken to be the $15$ closest points, including the center point. Further, $h$ is taken as the maximum distance between a center point $i$, and all its neighbouring points $j\in S_i$.

The point cloud is taken to be quasi-regularly spaced in this example only to determine the numerical order of convergence. Further, a small time step of $\Delta t = 0.1 h^2$ is used for the same reason. For a numerical solution $u$, relative errors in the solution are measured at $t=0.3$ as

\begin{equation}
	\epsilon_2 = \left[ \frac{\sum_{i=1}^N \left( u_i - u_{exact}(\vec{x}_i) \right)^2}{\sum_{i=1}^N \left( u_{exact}(\vec{x}_i) \right)^2 } \right]^{\frac{1}{2}}\,.
\end{equation}

A Crank--Nicolson time integration method is used. We consider two cases here: one with monomials up to the second order, and one with monomials up to the third order. The plots of the relative error against the number of points $N$ are shown in Figure~\ref{Fig:SurfaceHeatEquation}. In each of the cases, the experimental order of convergence matches the theoretical expectation. Due to the use of small time steps, the BiCGSTAB solver converged within very few iterations. Up to a tolerance of $10^{-10}$, the maximum number of iterations required for convergence was $5$, $5$, $4$ and $3$ for the second order case; and $3$, $3$, $3$ and $2$ for the third order case. In both cases, the number of iterations are reported in the order of increasing $N$.

To obtain an order of accuracy higher than third order, the size of the neighbourhoods would need to be increased. To use monomials up to order $p$, for a $k$ dimensional manifold, the number of monomial functions needed is given by the binomial coefficient ${p+k}\choose{k}$. Thus, to ensure that the local least squares systems for computing the differential operators are solvable, a necessary condition is that the number of neighbours should be more than ${p+k}\choose{k}$. The increased support sizes results in denser linear systems. This is one of the limitations of using classical GFDMs to approximate derivatives. To overcome this, modification in volumetric GFDMs to obtain higher order accuracy \cite{Milewski2012,Trask2017} would need to be extended to surfaces.

\begin{figure}
  \centering

  \includegraphics[width=0.7\textwidth]{./SurfaceHeatEqn_Sphere_N}

  \caption{Surface heat equation on a sphere: Relative errors against number of points $N$ in the domain. Second and third order convergence rates are marked with dashed lines. Numerical errors are shown by solid lines. The superscript of $\mathcal{P}^p$ indicates the maximum order of monomials being used in the computation of the differential operators.}

  \label{Fig:SurfaceHeatEquation}
\end{figure}

\subsection{Diffusion on a Torus with forcing}

We consider the diffusion equation again, but with forcing

\begin{equation}
	\label{Eq:HeatEquationF}
	u_t = \Delta_{M} u + f(u, \vec{x}, t)\,.
\end{equation}

The domain is taken to be a torus given by

\begin{equation}
	\left( 1 - \sqrt{ x^2 + y^2 } \right)^2 + z^2 = \frac{1}{9}\,.
\end{equation}

A manufactured solution is considered, as done in \cite{Fuselier2013}. The exact solution is given by

\begin{equation}
\label{Eq:TorusExact}
	u_{exact}(\vec{x}, t) = \frac{1}{8}\exp{(-5t)}x\left(x^4 - 10x^2y^2 + 5y^4 \right)\left( x^2 + y^2 - 60z^2 \right)\,.
\end{equation}

The forcing function is taken such that Eq.\,\eqref{Eq:TorusExact} satisfies Eq.\,\eqref{Eq:HeatEquationF}, and is given in \cite{Fuselier2013}. Once again, a Crank--Nicolson time integration method is used with a small $\Delta t$, as done above. From this example onwards, all points clouds are taken to be irregularly spaced, and are set up in a manner similar to \cite{Drumm2008}. The plots of relative error against the number of points $N$ are shown in Figure~\ref{Fig:SurfaceHeatEquation}. The figure illustrates that projection to the tangent plane along the central normal~(as done in this paper) is more accurate in the present context than the projection along the neighbour normals~(as done in some mesh-based methods \cite{Demanet2006}). Further, the experimental order of convergence is seen to match the theoretical expectation of second order for both types of projection. Up to a tolerance of $10^{-10}$, the maximum number of iterations required for convergence was $7$, $7$, $6$ and $5$, in order of increasing $N$. We note that the larger time steps used in the coming examples require higher number of iterations for convergence of the sparse linear system solver.

\begin{figure}
  \centering

  \includegraphics[width=0.7\textwidth]{./SurfaceHeatEqn_Torus_N}

  \caption{Diffusion on a torus with forcing: Relative errors against the number of points $N$ in the domain. Second order convergence rate is marked in blue. Numerical errors are shown in red~(projection to tangential plane along the central normal, as done in this paper), and in black~(projection along neighbour normals, as done in some mesh-based work \cite{Demanet2006}).}

  \label{Fig:SurfaceHeatEquationTorus}
\end{figure}

\subsection{Four Strip Problem on a Surface}

We now consider an elliptic problem with rough coefficients

\begin{equation}
\label{Eq:DivEtaGradPhi}
- \nabla_{M}\cdot\eta\nabla_{M} \phi = f\,,
\end{equation}

with appropriate boundary conditions. A common case to test numerical schemes is the volumetric equivalent of Eq.\,\eqref{Eq:DivEtaGradPhi} with discontinuous $\eta$ (for example, \cite{Trask2017}). Here, we consider not just rough coefficients, but $\eta$ with large jumps, up to several orders of magnitude.

Eq.\,\eqref{Eq:DivEtaGradPhi} with $f\equiv 0$ is solved on a wave surface defined by $z=\sin(2x)\sin(y)$, $x\in[0,4\pi]$, $y\in[0,4\pi]$, which is shown in Figure~\ref{Fig:DivEtaGradProblem_Domain}. The domain has four strips along the $x$ direction, each of the same width, with different values of the diffusion coefficient $\eta$. Starting from $x=0$, $\eta_1=10^4$, $\eta_2=10^0$, $\eta_3=10^2$, and $\eta_4 =10^0$. Figure~\ref{Fig:DivEtaGradProblem_Domain} shows the different strips for $\eta$. It also illustrates that the points are unevenly distributed, and that extra points are \textit{not} added along the interfaces where $\eta$ changes.

\begin{figure}
  \centering

  \includegraphics[width=0.7\textwidth]{./FourStripProblem_Domain}

  \caption{Discretized domain for the elliptic four-strip problem given by $z=\sin(2x)\sin(y)$. The colour indicates the different values of $\eta$.}

  \label{Fig:DivEtaGradProblem_Domain}
\end{figure}

Dirichlet boundary conditions are enforced at two ends, and pseduo-Neumann boundary conditions on the other two.

\begin{align}
\phi &= 0\phantom{abc} \text{if }\, x=0\,,\\
\phi &= 1\phantom{abc} \text{if }\, x=4\pi\,,\\
\vec{\hat{y}}\cdot\nabla_{M} \phi &= 0\phantom{abc} \text{if }\, y=0\,,\\
\vec{\hat{y}}\cdot\nabla_{M} \phi &= 0\phantom{abc} \text{if }\, y=4\pi\,,
\end{align}

where $\vec{\hat{y}}=(0,1,0)$. The domain is discretized with a total of $N=13918$ points. Due to the large jumps in $\eta$, the addition of the extra test functions for the computation of the diffusion operator, as explained in Section~\ref{sec:DiffusionLargeJumps}, is essential to prevent excessive numerical oscillations. Using these, the numerical solution is plotted in Figure~\ref{Fig:DivEtaGradProblem_Solution}, where the domain is warped by a scalar multiple of the solution $\phi$. The figure illustrates that the result shows a good agreement with the expectation that the solution does not vary in the $y$ direction, and that the linear slope of the solution in the $x$ direction is inversely related to the values of $\eta$.

\begin{figure}
  \centering

  \includegraphics[width=0.48\textwidth]{./Warp6x}
  \includegraphics[width=0.48\textwidth]{./Warp6x_OnlyHeight}

  \caption{Results for the elliptic four-strip problem: Domain warped by the solution. $x$ vs. $y$ vs. $z+6\phi$~(left) and $x$ vs. $y$ vs. $6\phi$~(right). The colour of the points indicates the value of $\phi$.}

  \label{Fig:DivEtaGradProblem_Solution}
\end{figure}

\subsection{Advection on a manifold}

Having considered a parabolic and an elliptic PDE, we now consider the hyperbolic problem of the transport equation

\begin{equation}
	\label{Eq:Advection}
	\frac{\partial \phi}{\partial t} + \vec{v}\cdot \nabla_{M}\phi = 0\,,\
\end{equation}

on the surface of a cone given by

\begin{equation}
 x^2 + y^2 = \frac{4}{9}z^2\,,
\end{equation}

with $-6 \leq z \leq 0$. To discretize Eq.\,\eqref{Eq:Advection}, we once again take advantage of the fact that developments in volumetric GFDMs can be directly carried over here. Numerical methods to discretize the advection term for the volumetric case have been widely studied, including using meshfree GFDMs \cite{Praveen2007, SeifarthThesis}. Here, we follow the work of Seifarth \cite{SeifarthThesis} for volumetric meshfree GFDMs. A upwind discretization of Eq.\,\eqref{Eq:Advection} leads to the following semi-discrete form

\begin{equation}
	\frac{d\phi_i}{dt} = -2 \sum_{\substack{j\in S_i \\ j \neq i}} c_{ij}^{\vec{v}\cdot\nabla_M} (\phi_{ij} - \phi_i) \,,
\end{equation}

where

\begin{align}
	c_{ij}^{\vec{v}\cdot\nabla_M}  &= c_{ij}^{M,x}v_i^1 + c_{ij}^{M,y}v_i^2 + c_{ij}^{M,z}v_i^3\,,\\
	\phi_{ij} &= \frac{1}{2} \left[ \left( 1 + \text{sign}(\delta\vec{x}_{ij} \cdot \vec{v}_i ) \right)\phi_{ij}^{+}  +  \left( 1 - \text{sign}(\delta\vec{x}_{ij} \cdot \vec{v}_i ) \right)\phi_{ij}^{-} \right]\,,
\end{align}

where $\delta\vec{x}_{ij} = \vec{x}_j - \vec{x}_i$ are the distances computed in $\mathbb{{R}}^3$; and $\phi_{ij}^{+}$ and $\phi_{ij}^{-}$ are the reconstructed values. Here, we use a MUSCL reconstruction  with a Superbee limiter. A pure upwind scheme is also shown for comparison. For both cases, a SDIRK2 \cite{Alexander1977, Ellsiepen1999} implicit second order method of time integration is used, which enables the use of large time steps. A linearization of the arising sparse implicit system is done before using a BiCGSTAB solver. More details of the scheme, including the time integration, can be found in \cite[Section 4.3]{SeifarthThesis}.

A Gaussian bell is transported on the cone. The initial condition is given by

\begin{equation}
	\label{Eq:Advection_IC}
	\phi(\vec{x},0) = \left\{ \begin{array}{ll}
         \frac{ \exp\left( - \| \vec{x} - \vec{x}_0\|^2\right) - \exp\left(-25\right) }{1 - \exp\left( -25 \right) } &\; \text{if } \| \vec{x} - \vec{x}_0\|^2 < 5^2\,, \\
          0,      &\; \text{elsewhere }\,,
        \end{array} \right.
\end{equation}

for $ \vec{x}_0 = (2,0,-3)$. The velocity field is taken to be $\vec{v} = (-y,x,0)$. Note that this is not divergence-free on the manifold, and that it lies on the tangent bundle of the manifold.

Time integration is performed with a time step of $\Delta t = 0.02$ until $t_{end} = 2\pi$, which corresponds to one full rotation. The domain is discretized with $N=52\,193$ points. The results using the scheme mentioned above are shown in Figure~\ref{Fig:Advection_Results}. Results for a pure upwind scheme~(without any reconstruction) with the same time integration method and same initial conditions are shown in Figure~\ref{Fig:Advection_PureUpwind}. These figures illustrate that the MUSCL reconstruction with Superbee limiter results in minimal numerical diffusion, while the pure upwind scheme causes excessive numerical diffusion.
Further, a direct central difference approximation of the $\vec{v}\cdot\nabla$ operator leads to unstable simulations. Each of these observations agrees with the expectation for volumetric advection. The errors can be quantified as follows. For the numerical domain with $N=52\,193$, as used in Figures~\ref{Fig:Advection_Results} and \ref{Fig:Advection_PureUpwind}, at $t=0$, $\text{max}(\phi) = 0.99$ for both simulations (the maximum value of $1$ of Eq.\,\eqref{Eq:Advection_IC} is not attained as no point is present exactly at $\vec{x}_0$). At $t=2\pi$, $\text{max}(\phi) = 0.36$  for the pure upwind case, and $\text{max}(\phi) = 0.97$ for the MUSCL and Superbee case.

\begin{figure}
  \centering

  \includegraphics[width=0.3\textwidth]{./ConeAdvection0000}\phantom{abcdefghij}
  \includegraphics[width=0.3\textwidth]{./ConeAdvection0011}\\
  \includegraphics[width=0.3\textwidth]{./ConeAdvection0022}\phantom{abcdefghij}
  \includegraphics[width=0.3\textwidth]{./ConeAdvection0033}\\

  \caption{Advection on the surface of a cone for $N=52\,193$: evolution of the solution $\phi$ at different times using a MUSCL reconstruction with a Superbee limiter. $t=0 s$~(top left), $t=0.66 s$~(top right), $t=1.32 s$~(bottom left), and $t=1.98 s$~(bottom right).}

  \label{Fig:Advection_Results}
\end{figure}

\begin{figure}
  \centering

  \includegraphics[width=0.3\textwidth]{./ConeAdvection_PureUpwind0022}\phantom{abcdefghij}
  \includegraphics[width=0.3\textwidth]{./ConeAdvection_PureUpwind0033}

  \caption{Advection on the surface of a cone for $N=52\,193$: evolution of the solution $\phi$ at different times for a pure upwind scheme. $t=1.32 s$~(left), and $t=1.98 s$~(right).}

  \label{Fig:Advection_PureUpwind}
\end{figure}

After, one full rotation, at $t=2\pi$, the analytical solution matches the initial condition. We measure relative errors for this with varying $N$. The plots of the same are shown in Figure~\ref{Fig:AdvectionErrors}. A smaller time step is used in these simulations $\Delta t = 0.1h$. The smoothing lengths considered in the simulations for Figure~\ref{Fig:AdvectionErrors} start with $h=0.6$ and are consecutively halved till $h=0.075$, which results in the number of points varying from $N=2960$ to $N=853837$. The value of the peak of the Gaussian is also shown in the same figure to illustrate the numerical diffusion of the peak of the solution.

\begin{figure}
  \centering

  \includegraphics[width=0.47\textwidth]{./AdvectionErrors1}
  \includegraphics[width=0.47\textwidth]{./AdvectionErrors2}

  \caption{Advection on the surface of a cone: Relative errors~(left) and the error in the peak of the solutions~(right) at $t=2\pi$ for varying $N$. For the pure upwind case~(black), and the MUSCL + Superbee case~(red).}

  \label{Fig:AdvectionErrors}
\end{figure}

This example, once again, illustrates the ease of taking over volume based GFDM ideas to the surface PDE discretization setting explained in this paper. In fact, for the simulations in this section, not only was the work of volumetric methods carried over, the code of Seifarth \cite{SeifarthThesis} was also used directly.

\subsection{Cahn--Hilliard Equation}

We now consider the PDE system of the Cahn-Hilliard equation \cite{Cahn1958} on a closed manifold (with no boundaries). It describes the process of phase separation \cite{Kim2016}. Numerical solutions to the equation on manifolds has had a lot of interest \cite{Chen2014, Du2011}. Here, we consider the Cahn-Hilliard equation as a pair of coupled second order PDEs, as done in \cite{Gera2017}.

\begin{align}
	\frac{\partial f}{\partial t} &= \frac{1}{Pe} \Delta_{M} \mu \label{Eq:CahnHilliard1}\,,\\
	\mu &= \frac{dg}{df} - Cn^2 \Delta_{M}f\label{Eq:CahnHilliard2}\,,
\end{align}

where $f$ is an order parameter on the manifold, $\mu$ is the chemical potential, $Pe$ is the surface Peclet number, and $Cn$ is the Cahn number. As commonly done, the function $g$ is taken according to a double well potential $g(f) = \frac{f^4}{4} - \frac{f^2}{2}$.

Eqs.~\eqref{Eq:CahnHilliard1} and \eqref{Eq:CahnHilliard2} are solved on a surface of an airplane, as shown in Figure~\ref{Fig:Airplane_Domain}. The airplane CAD files are based on that by NASA CRM \cite{Vassberg2008}. Note that each part of the airplane, including the wings, have non-zero thickness, and there is a layer of points on each side of the geometry. Thus, the discretized airplane considered here has no boundaries, and forms a closed manifold.

\begin{figure}
  \centering

  \includegraphics[width=0.5\textwidth]{./AirplaneDomain}

  \caption{Discretized domain for the Cahn-Hilliard equation.}

  \label{Fig:Airplane_Domain}
\end{figure}

This example also illustrates that the present method can handle complex geometries and non-smooth surfaces. The connection between the wings and the fuselage, for instance, has sharp changes. A first order time-integration scheme is used

\begin{align}
	\frac{f^{(n+1)} - f^{(n)}}{\Delta t} &= \frac{1}{Pe} \Delta_{M} \mu^{(n+1)}\,,\\
	\mu^{(n+1)} &= \frac{dg}{df}(f^{(n)}) - Cn^2 \Delta_{M}f^{(n+1)}\,,
\end{align}

which is solved in one large coupled linear system. The evolution of $f$ is shown in Figure~\ref{Fig:CahnHilliard_Results} for $Pe = 1$, $Cn = 0.5$, $\Delta t = 10^{-4}$and $N= 15\,713$ points in the domain. Our simulation results show similar trends to those exhibited in mesh-based solutions \cite{Chen2014,Du2011}. The concentration $f$ first evolves to the range $[-1,1]$, and the separation continues until a steady state is reached.

\begin{figure}
  \centering

  \includegraphics[width=0.48\textwidth]{./Airplane_CahnHilliard_0}
  \includegraphics[width=0.48\textwidth]{./Airplane_CahnHilliard_15}\\
  \includegraphics[width=0.48\textwidth]{./Airplane_CahnHilliard_50}
  \includegraphics[width=0.48\textwidth]{./Airplane_CahnHilliard_225}\\
  \includegraphics[width=0.48\textwidth]{./Airplane_CahnHilliard_375}
  \includegraphics[width=0.48\textwidth]{./Airplane_CahnHilliard_500}\\

  \caption{Cahn-Hilliard Equation on the surface of an airplane: evolution of the order parameter $f$ at different times.}

  \label{Fig:CahnHilliard_Results}
\end{figure}

\section{Conclusion}
\label{sec:Conclusion}

We presented a new meshfree approach to solving PDEs defined on manifolds embedded in $\mathbb{R}^n$. This approach is applicable for manifolds both with and without boundaries. Differential operators are computed directly on the tangent space, and they can be applied directly to function values on the manifold. The discretized domain consists of points only on the manifold, with no bulk discretization around the surface. Consequently, the method scales with the actual dimension of the manifold, and \textit{not} the dimension of the embedding space.

One of the biggest advantage of this method lies in the fact that the differential operators are computed on the tangent space in the same way as regular, volume-based, meshfree methods for the dimension of the tangent space. As a result, all developements in classical volumetric meshfree GFDMs can be directly carried over to manifolds. This was illustrated with different examples, for the optimization of the Laplacian stencil, for treating diffusion with large jumps in the diffusion coefficients, and for discretizing advection terms.

The applicability of this method was shown with numerical examples for a parabolic PDE, an elliptic PDE, and a hyperbolic PDE. Numerical results showed good compatibility with analytical solutions. The numerical simulations also showed that this method can handle unevenly distributed particles, and sharp edges in the point cloud. The simulations showed promising results, and suggest that this method could be extended to solve a wider range of problems on surfaces.

\section*{Ackowledgements}

The first author would like to thank Dr. Tobias Seifarth for his input on advection methods in meshfree GFDMs, and for sharing his code on the same.

\appendix
\section{Optimized Surface Laplacian Stencils}
\label{App:Laplace}

The Laplace stencil is split into two parts. One which satisfies the consistency conditions of Eq.\,\eqref{Eq:TPL_Consistency}, and one which pushes the stencil towards positivity.

\begin{equation}
\label{Eq:LaplaceSplit}
	c_{ij_{T}}^{\Delta_T} = \sigma_{ij_{T}}^{\Delta_{T}} + \alpha^{\Delta_{T}}d_{ij_{T}}^{\Delta_{T}}\,,
\end{equation}

where the $\sigma_{ij_{T}}^{\Delta_{T}}$ satisfy the consistency conditions, and $d_{ij_{T}}^{\Delta_{T}}$ are used to improve the conditioning of the Laplace stencil. They are given by

\begin{align}
	\sum_{j\in S_i}\sigma_{ij_{T}}^{\Delta_{T}}m_{j_{T}} &= \Delta_{T} m(\vec{x}_i) \qquad \forall m\in\mathcal{P}_{T}\,,\\
	\sigma_{ii}^{\Delta_{T}} &= \tilde{A}_c \,,\\
	\text{min } J_i &= \sum_{j\in S_i} \left( \frac{\sigma_{ij_{T}}^{\Delta_{T}}}{W_{ij_{T}}} \right)^2\,,
\end{align}

where $\tilde{A}_c \notin \{0,1\}$ is some fixed central stencil value. And

\begin{align}
	\sum_{j\in S_i}d_{ij_{T}}^{\Delta_{T}} m_{j_{T}} &= 0 \qquad \forall m\in\mathcal{P}_{T}\,,\label{Eq:DD_dStencil_1}\\
	d_{ii}^{\Delta_{T}} &= 1 \,,\label{Eq:DD_dStencil_2}\\
	\text{min } J_i &= \sum_{j\in S_i} \left( \frac{d_{ij_{T}}^{\Delta_{T}}}{W_{ij_{T}}}  \right)^2\,.\label{Eq:DD_dStencil_3}
\end{align}

We note that the stencil coefficients $\sigma$ and $d$ can be computed by one minimization with different right hand side vectors \cite[Section A.4]{Suchde2018_Thesis}, and thus, performing this procedure does not increase the computation time significantly.  Now, $\alpha^{\Delta_{T}}$ is computed by minimizing the functional

\begin{equation}
	g^{\Delta_{T}} = \frac{\sum_{j \in S_i} \left(c_{ij_{T}}^{\Delta_{T}} \right)^2}{\left (c_{ii_{T}}^{\Delta_{T}} \right)^2} \,.
\end{equation}

Setting $\frac{\partial g^{\Delta_{T}}}{\partial c_{ii}^{\Delta_{T}}} = 0$, we get

\begin{equation}
	\label{Eq:DD_alpha}
	\alpha^{\Delta_{T}} = \frac{\langle\vec{\sigma},\vec{d} \;\rangle \tilde{A}_c - \langle\vec{\sigma},\vec{\sigma}\;\rangle }{ \langle\vec{\sigma},\vec{d}\;\rangle - \langle\vec{d},\vec{d}\;\rangle \tilde{A}_c}\,,
\end{equation}

where
\begin{equation}
	\langle \vec{\sigma},\vec{d} \;\rangle = \sum_{j \in S_i} \sigma_{ij_{T}}^{\Delta_{T}} d_{ij_{T}}^{\Delta_{T}}\,.
\end{equation}
The proof of the above has been shown in \cite{Suchde2018_Thesis} (Section 2.5.5 and Appendix A.4). Using Eq.\,\eqref{Eq:DD_alpha}, we can compute the numerical tangent plane Laplacian according to Eq.\,\eqref{Eq:LaplaceSplit}. This, in turn, gives the numerical surface Laplacian.

\section{Extra Consistency Conditions on the Anisotropic Laplacian Operator}
\label{App:Diffusion}

Here, we prove that the extra consistency conditions Eq.\,\eqref{Eq:DiffExtra1} -- Eq.\,\eqref{Eq:DiffExtra3} imposed on the diffusion operator to deal with large jumps in the diffusion coefficient are valid. To simplify notation, we show the same for the volumetric $2$ dimensional case. And, without loss of generality,  we assume that $\vec{s}$ is along the $x$ direction. They are derived by some algebraic manipulations to improve numerical conditioning. Further, these procedures also avoid the numerical computation of derivatives of $\kappa$, and replace it by the computation of derivatives of $\log(\kappa)$.

Firstly, Eq.\,\eqref{Eq:DiffExtra1} is obtained by the addition of the test function $\frac{1}{\kappa}$, the action of the diffusion operator on which can be given by

\begin{align}
\nabla \cdot \left( \kappa\nabla \left( \frac{1}{\kappa} \right) \right) &= \nabla\cdot \left( - \frac{\nabla \kappa}{\kappa} \right) \,,\\
&= -\nabla \cdot \left( \nabla \left(  \log\kappa     \right)\right)\,, \\
&= -\Delta \left( \log\kappa \right)\,.
\end{align}

Eq.\,\eqref{Eq:DiffExtra2} is obtained by the addition of the test function $\frac{\delta x}{\kappa}$, the action of the diffusion operator on which can be given by

\begin{align}
\nabla \cdot \left( \kappa\nabla \left( \frac{\delta x}{\kappa} \right) \right) &= \nabla\cdot \left( \nabla\left( \delta x \right) - \delta x\frac{\nabla \kappa}{\kappa} \right) \,,\\
&= \Delta \left( \delta x \right) - \nabla\left( \delta x \right)\cdot \nabla \left(\log\kappa \right) - \delta x \Delta \left( \log\kappa \right)\,\\
&= 0 - \frac{\partial}{\partial x} \left( \log\kappa \right) - 0\,,
\end{align}

noting that $\delta x$ evaluated at the central point $i$ is $0$. Lastly, Eq.\,\eqref{Eq:DiffExtra3} is obtained by the addition of the test function $\frac{\delta x ^2}{\kappa}$, the action of the diffusion operator on which can be given by

\begin{align}
\nabla \cdot \left( \kappa\nabla \left( \frac{\delta x^2}{\kappa} \right) \right) &= \nabla\cdot \left( \nabla\left( \delta x \right)^2 - \delta x^2\frac{\nabla \kappa}{\kappa} \right) \,,\\
&= \Delta \left( \delta x \right)^2 - \nabla\left( \delta x \right)^2\cdot \nabla \left(\log\kappa \right) - \delta x^2 \Delta \left( \log\kappa \right)\,\\
&= 2 -0 - 0\,.
\end{align}