\section{Introduction}

Boris (Boaz) Abramovich Trakhtenbrot (1921-2016) was one of the pioneers
in recognizing the usefulness of Monadic Second Order Logic $\MSOL$
for treating situations in automata theory \cite{trakhtenbrot1962finite}.
In \cite{ar:FischerMakowsky08}, published in the Festschrift for Boaz' 85th birthday,
Eldar Fischer and the second author gave an application of Monadic Second Order Logic
to graph polynomials, by proving that for a  graph polynomial $P(G;\bar{x})$ definable in $\MSOL$
and a sequence of recursively defined graphs $G_i: i \in \N$
the sequence of polynomials $P(G_i;\bar{x})$ satisfies a linear recurrence relation over the
polynomial ring $\Z[\bar{x}]$.

Ernst Specker (born in 1920-2011), together with Christian Blatter (1935-2021), was the first to use
Monadic Second Order Logic to prove a meta-theorem in counting combinatorics, \cite{pr:BlatterSpecker84}.
In order to celebrate Boaz' centenary we shall give another application of Monadic Second Order Logic
to $\MSOL$-definable  graph polynomials $P(G;\bar{x})$, where the proof again uses the Specker-Blatter Theorem.
We fix $\bar{a} \in \N^r$ and look at the sequence  $p(n, \bar{a})=P(K_n;\bar{a})$ modulo $\mu \in \N$,
where $K_n$ is the complete
graph on $n$ vertices. We show that, under simple assumption on $\bar{a}$ and $\mu$,
this sequence is ultimately periodic modulo $\mu$.

Our results
are easy, but possibly unexpected,
applications of the Specker-Blatter Theorem.
They illustrate more the power of it,
by using general methods of Monadic Second Order Logic $\MSOL$, rather than
applying combinatorial arguments specially tailored to a particular case.

\subsection{Some graph polynomials}

Let $G =(V(G), E(G))$ be a finite graph.
We put $n(G) = |V(G)|$, $m(G) = |E(G)|$, $\kappa(G)$ is the number of connected components of $G$.
From a logical point of view $G$ can be represented in various ways.
The vocabulary $\tau_{graph}$ consists of one binary relation symbol for $E(G)$
and the vocabulary $\tau_{hgraph}$ consists of two unary predicates, one for vertices $V(G)$ and one for edges $E(G)$
and  one binary relation symbol for $R(G)$ for the incidence relation between edges and vertices.
$\tau_{hgraph}$ is also suitable for representing hypergraphs.
{\em Induced subgraphs} are substructures in the case of $\tau_{graph}$.
and {\em subgraphs} are substructures in the case of $\tau_{hgraph}$.
In other words, let $A \subseteq V(G)$ be a set of vertices.
The {\em induced subgraph} of $G$
generated by $A$ is the graph
$G[A] = (A, E(G) \cap A^2)$.
On the other hand the graph $(A,E)$ is a subgraph of $G$ for any
$F \subseteq E(G)$.

$\MSOL$ on graphs allows quantification over subsets of vertices only.
$\MSOL$ on hypergraphs allows quantification over subsets of vertices and subsets of edges.
$\MSOL$ on hypergraphs has the same expressive power as Second Order Logic $\SOL$ on graphs where second order quantification
is {\em restricted to unary predicates and binary relations which are subsets of the edge relation}.
We denote this version by $\GMSOL$, for Guarded Monadic Second Order Logic.
We assume the reader is familiar with Monadic Second Order Logic $\MSOL$.
For the Monadic Second Order theory of graphs the reader is referred to the encyclopedic \cite{bk:CourcelleEngelfriet2011}.
$\CMSOL$ is the logic extending $\MSOL$ with modular counting quantifiers.

A {\em graph polynomial $P(G,\bar{X})$} is a {\em function $P$} which associates with
a graph $G$ a polynomial $P(G,\bar{X}) \in \Z[\bar{X}]$, and which is {\em invariant}
under graph isomorphisms (disregarding the labels).
Here $\bar{X} =(X_1, \ldots , X_s)$.

\medskip
Let us look at some examples:
\begin{enumerate}[(i)]
\itemsep=0.9pt
\item
Let $i_k(G)$ denote the number of independents sets $A \subseteq V(G)$ of size $k$.
\\
The {\em independence polynomial $In(G,X)$} is defined as
$$
In(G,X)= \sum_{k=0}^{n(G)} i_k(G) X^k
$$
\item
Let $c_k(G)$ denote the number of sets $A \subseteq V(G)$ which induce a clique of size $k$.
\\
The {\em clique polynomial $Cl(G,X)$} is defined as
$$
Cl(G;X) =\sum_{k=0}^{n(G)} c_k(G) X^k
$$
\item
Let $\lambda \in \N$. $\chi(G,\lambda)$ denotes the number of proper $\lambda$-colorings of $G$.
\\
By the well-known observation of G. Birkhoff (1912), $\chi(G,\lambda)$ is a polynomial in $\N[\lambda]$.
\\
We denote by $\chi(G,X)$ the extensions of $\chi(G,\lambda)$ to a polynomial over the complex numbers,
$\chi(G, X) \in \C[X]$.
$\chi(G,X)$ is called the {\em chromatic polynomial of $G$}.
\item
The {\em Tutte polynomial $T(G;X,Y)$} is defined as
\begin{gather}
T(G;X,Y) = \sum_{A \subset E(G)} (X-1)^{\kappa(A) - \kappa(G)} \cdot (Y-1)^{|A| + \kappa(A) -|V(G)|}
\label{tutte-def}
\end{gather}
where $\kappa(S)$ is the number of connected components of the spanning subgraph $G[S]=(V(G), S)$.
\item
The {\em matching polynomials} come in two versions:
Let $G$ be a graph and $m_k(G)$  be the number of matchings of size $k$ of $G$.
The {\em generating matching polynomial} $M(G:X)$ of $G$ is defined as
$$
M(G;X) = \sum_{k=0}^{\lfloor n/2 \rfloor} m_k(G) X^k
$$
and the {\em matching defect} aka {\em acyclic polynomial} $\alpha(G;X)$ is defined as
$$
\alpha(G;X) = \sum_{k=0}^{\lfloor n/2 \rfloor} (-1)^k m_k(G) X^{n-2k}
$$
The two are related by the equations
$$
\alpha(G;X) = X^n M(G;-X^{-2}) \text{  and  }
M(G;X)=  (-i)^n X^{n/2} \alpha(G;i X^{-1/2})
$$
There is also a bivariate version
$$
\bar{M}(G;X,Y) = \sum_{k=0}^{\lfloor n/2 \rfloor} (X)^k m_k(G) Y^{n-2k}
$$
where $\bar{M}(G;-1,Y)=\alpha(G;Y)$ and $\bar{M}(G;X,1)=M(G;X)$.
\item
In \cite{averbouch2008most,averbouch2010extension}
the authors introduce a most general edge elimination polynomial in three indeterminates $\xi(G;X,Y,Z)$.
\begin{gather}
\xi(G;X,Y,Z) = \sum_{(A \sqcup B) \subseteq E}
X^{\kappa(A \sqcup B) - c(B)}
\cdot
Y^{|A|+|B| - c(B)}
\cdot
Z^{c(B)}
\label{xi-def}
\end{gather}
Here $c((V,E))$ is the number of connected components of $(V,E)$
which have at least one edge.
Both the matching polynomials and the Tutte polynomial are substitution instances of $\xi(G;X,Y,Z)$.
In \cite{trinks2011covered,trinks2012proving}
other trivariate graph polynomials are discussed which are equivalent to $\xi(G;X,Y,Z)$:
Among them the subgraph counting polynomial $S(G,X,Y,Z)$ and the covered components polynomial $C(G;X,Y,Z)$.
All these trivariate graph polynomials are
substitution instances of each other.
The subgraph counting polynomial $S(G;X,Y,Z)$
is defined as\smallskip
$$
S(G;X,Y,Z) = \sum_{H=(W,F) \subseteq G}
X^{|W|}
Y^{\kappa(H)}
Z^{|F|}
$$
where $H$ ranges over all subgraphs of $G$.

\vspace{1mm}
The covered components polynomial $C(G;X,Y,Z)$ is defined as
$$
C(G;X,Y,Z) = \sum_{A\subseteq E}
X^{\kappa(([n],A))}
Y^{|A|}
Z^{c(([n],A))}
$$
\end{enumerate}
In \cite{trinks2011covered} it is shown how  $\xi(G;X,Y,Z)$ and $C(G;X,Y,Z)$ are related:
\begin{proposition}
\label{trinks}
$$
C(G;X,Y,Z) = \xi((G;X,Y, XYZ-XY)
\text{   and   }
\xi((G;X,Y,Z) = C(G;X,Y, \frac{Z}{XY} +1)
$$
\end{proposition}
A similar relation is given for $\xi(G;X,Y,Z)$ and $S(G;X,Y,Z)$ in \cite{trinks2012proving} which is a bit more complicated,
but no needed for this paper.

\begin{remark}
We note that both the Tutte polynomial and $\xi(G;,X,Y,Z)$ involve negative exponents,
whereas $S(G;X,Y,Z)$ and $C(G;X,Y,Z)$ do not.
This is the source of the difference between Theorem \ref{th:1} and Theorem \ref{th:4}.
\label{negative-exp}
\end{remark}

\subsection{The case of $G =K_n$}

Here we shall be concerned with computing a graph polynomial $P(G;\bar{X})$ in $k$ indeterminates
for the case where the graph $G$ is $K_n$, the complete graph on $n$ vertices.
We define, for
fixed non-negative integers $\bar{b} \in \N^k$
the sequence
$$P_{n}(\bar{b})= P(K_n,\bar{b})$$
for fixed values $\bar{b} \in \Z^k$.
Can we make some general statement about the sequence $P_n(\bar{b})$?
In some cases computing $P_n(\bar{b})$ is very easy, using trivial observations or simple recurrence relations.
However, the resulting graph polynomials may be unexpectedly complicated.
The Tutte polynomial and the
matching polynomial will illustrate this in the sequel.

\medskip
We compute $P_n(\bar{b})$ first for some straightforward cases:
\begin{enumerate}[(i)]
\itemsep=0.85pt
\item
{\em Independence polynomial}:
$$
i(K_n,k) =
\begin{cases}
0 & k \geq 2 \\
n & k=1 \\
1 & k=0
\end{cases}
$$
\eject

\noindent Hence,
$$
In(K_n, b) = i_1(K_n) \cdot b + c_0(K_n) = nb +1
$$
\item
{\em Clique polynomial}: $c_k(K_n) = {n \choose k}$, hence
$$
Cl(K_n,b) = \sum_{k=0}^n {n \choose k} b^k = (b+1)^n
$$
\item
{\em Chromatic polynomial}:
$\chi(K_n,b) = 0$ for $n\geq b+1$.
\end{enumerate}

\subsection{The Tutte polynomial}

The case of the Tutte polynomial is a bit more complicated, but of special interest.
We note that
\begin{enumerate}[(i)]
\itemsep=0.85pt
\item
$T(K_n;2,1)$
counts the number of forests on $n$ vertices.
\item
$T(K_n;1,1)$
counts the number of trees on $n$ vertices.
\item
$T(K_n;1,2)$
counts the number of connected graphs on $n$ vertices.
\end{enumerate}

The earliest computation of  $T_n(X,Y) = T(K_n,X,Y)$
can be found in \cite{tutte1954contribution} which already gives a recursive computation.
I. Gessel
\cite{gessel1979depth,gessel1995enumerative,gessel1996tutte}
and independently I. Pak \cite{pak-tutte}, proved:
$$
T(K_n; a,b) = \sum_{k=1}^n {(n-1) \choose (k-1)} \left(a + \sum_{i=1}^{k-1} b^i \right) T(K_{k-1};1,b) \cdot T(K_{n-k};a,b)
$$
I. Pak  in \cite{pak-tutte} also lists many other evaluations of $T(K_n; X,Y)$
with their combinatorial interpretations.

However,
in \cite{ar:BiggsDamerellSand72} it is noted that $T_n(X,Y)$
does not satisfy a linear recurrence which is independent of $n$.

\begin{proposition}[N.L. Biggs, R.M. Damerell and D.A. Sand]
There is no linear recurrence relation which computes the sequence $T_n(a,b)$ for fixed $a,b \in \Z$.
\end{proposition}

\subsection{The matching polynomial}

From \cite{godsil1981hermite} we know that for the defect matching polynomial (aka the acyclic polynomial)
$$
\alpha(K_n;X) = He_n(X)
$$
where $He_n(X)$ are the probabilist's Hermite polynomials for $n \in \N$.
The proof of this is due to C. Heilmann and E. Lieb
\cite{HeilmannLieb}.
From \cite[Equations 3.4 and 3.8]{carlitz1953congruence} one can derive\footnote{
Thanks to V. Rakita for checking this.
} that
the polynomials $He_n(X)$ satisfy the modular recurrence relation
$$
He_{n+m}(X) = He_n(X) \cdot  He_{m}(X) = He_n(X) \cdot X^m \left(\hspace*{-3mm}\mod{\mu}\right).
$$
and with $He_{0}(X) = 1$ and one gets  $He_m(X) =X^m$.

\begin{proposition}[Carlitz, 1953]
\label{pr:carlitz}
For every $a \in \N$ the sequence
$\alpha(K_n;a) = He_n(a) =a^n$ is ultimately period modulo $\mu$.
\end{proposition}
Our Theorem \ref{th:3} below shows the ultimate periodicity modulo $\mu$ for $M(K_n,a)$
without using the connection to the Hermite polynomials.

\subsection{Computing $P_n(\bar{b})$ modulo an integer $\mu$}

Let $\mu \in \N$.
We compute $P_n(\bar{b})$ modulo $\mu$ and observe:
\begin{enumerate}[(i)]
\itemsep=0.9pt
\item
For every $b \in \Z$ and $\mu \in \N$ the sequence
$In(K_n, b)= nb+1$ is {\em ultimately periodic}.
\item
For every $b \in \Z$ and $\mu \in \N$ the sequence
$Cl(K_n,b) =(b+1)^n$ is {\em ultimately periodic}.
\item
The sequence $\chi(K_n,b)$ is ultimately constant, hence  {\em ultimately periodic}.
\end{enumerate}
We shall see that this is the case
for a very large class of graph polynomials
subject to a definability condition in $\MSOL$.

Our main results are for the Tutte polynomial $T(G;X,Y)$, the bivariate matching polynomial $\bar{M}(G;X)$
and the trivariate edge elimination polynomial $\xi(G;X,Y,Z)$.
\begin{theorem}
\label{th:1}
For every $a,b, \mu \in\N^{+}$ with $a,b>1$,
$\gcd(a-1, \mu )=1$ and $\gcd(b-1, \mu )=1$, \\
the sequence $T(K_{n},a,b)$ is ultimately periodic modulo $\mu$.
\end{theorem}
Similarly, we also get for
the edge elimination polynomial $\xi(G;X,Y,Z)$,
the subgraph counting polynomial $S(G;X,Y,Z)$, and
the covered components polynomial $C(G;X,Y,Z)$:
\begin{theorem}
\label{th:4}
\begin{enumerate}[(i)]
\item
For every $a,b, c, \mu \in\N^+$
the sequences
$S(K_{n},a,b,c)$ and
$C(K_{n},a,b,c)$
are ultimately periodic modulo $\mu$ for every $\mu$.
\item
Assume  $ab$ divides $c$.
\\
Then $\xi(K_{n},a,b,c)$ is ultimately periodic modulo $\mu$ for every $\mu$.
\end{enumerate}
\end{theorem}
\begin{remark}
For $\xi(K_{n},a,b,c)$ one needs an additional condition like in Theorem \ref{th:1} due to the negative
exponents in the definition
\ref{xi-def}.
of $\xi(G;,X,Y,Z)$.
\end{remark}
\begin{theorem}
\label{th:3}
For every $a,b, \mu \in\N^{+}$
the sequence $\bar{M}(K_{n},a,b)$ is ultimately periodic modulo $\mu$.
\end{theorem}
Theorem \ref{th:3} can be generalized:
\begin{theorem}
\label{th:2}
If $P$ is a graph polynomial definable in $\GMSOL_{hgraph}$ without order,
then for every $\mu \in \N$ and every $a,b \in \Z$
the sequence $P(K_n, a,b)$ is ultimately periodic modulo $\mu$.
\end{theorem}
\begin{remark}
We note that Theorem \ref{th:1} is not a special case of Theorem \ref{th:2},
as the Tutte polynomial seems not to be $\GMSOL_{hgraph}$-definable without an order on the vertices.
\end{remark}

The three theorems only assert that the sequences are ultimately periodic modulo $\mu$, without
any indication of the length of the periodicity or the initial segment before the periodicity starts.

\subsection{Methods}

The proofs use tools from logic and combinatorics.
In particular
\begin{itemize}
\itemsep=0.9pt
\item
{\em Definability in Monadic Second Order Logic $\MSOL$};
\item
{\em Definability} in the extension $\CMSOL$ of $\MSOL$, where modular counting quantifiers are added;
\item
The {\em Specker-Blatter Theorem},
which gives a sufficient condition on when certain $\CMSOL$-definable density functions
are ultimately periodic if considered modulo $\mu$.
\end{itemize}

We assume that our readers are familiar with $\MSOL$ and first explain the {\em Specker-Blatter Theorem}.

\subsection{Related results}

The particular sequence
$T_n(1,b)=T(K_n;1,b)$
has been analyzed by
A. P. Mani and R.J. Stones,
\cite{ar:ManiStones2016},
for
$\mu = p^k$ where $p$ is an odd prime and $k \in \N$.

Let $\phi(n)$ be the {\em Euler totient function}
which is defined as the number of integers $a \in \{1, 2, . . . , m\} =[m]$ such that
$gcd(a, m) = 1$,  cf. \cite{bk:GrahamKnuthPatashnik94}.

\begin{proposition}[A.P. Mani and R.J. Stones]
Let $p$ be a prime, and let $k$ be a positive integer. For $b, n \in \Z$
such that $n \geq p^k$ and $b \neq 1 \mod{p}$, it holds that
$$
T_n(1,b) =\mod{p^k}
\begin{cases}
b^{\frac{\phi(p^k )}{2}} C_{n - \phi(p^k}(b) & \text{ if }  p \geq 3 \text{ and } n > p \\
b^{\frac{\phi(p)}{2}} - 1                   & \text{ if }  p \geq 3 \text{ and } n = p \\
1                                  & \text{ if }  p = n = 2\\
2                                  & \text{ if }  p = k = 2 \text{ and } n = 4\\
0                                  & \text{ otherwise}
\end{cases}
$$
\end{proposition}
This is much more informative than our Theorem \ref{th:1} for the case $a=1$ and $\mu = p^k$.
In \cite{ar:ManiStones2016} they formulate also a conjecture for $T(K_n;a,b)$ for general $a$, but still for $\mu = p^k$.

\begin{conjecture}[A.P. Mani and R.J. Stones]
Let $p$ be an odd prime and let $k$ be a positive integer.
\begin{enumerate}[(i)]
\itemsep=0.9pt
\item
If $n \geq p^k$, $a,b \in \Z$, and $b \not\equiv 1 \mod{p}$, then
modulo $p$ we have
$$
T(K_n;a,b) \equiv\mod{p}
\begin{cases}
b^{\frac{\phi(p)}{2} - 1} & \mbox{   if   } p\geq3, n=p  \\
& \mbox{ and } a=1 \mod{p}\\
b^{\frac{\phi(p^k)}{2}} \cdot T(K_{n - \phi(p^k)};a,b) & \mbox{otherwise}
\end{cases}
$$
\item
If $n \geq p^k$, $a,b \in \Z$, and $b \equiv 1 \mod{p}$, then
$$
T(K_n;a,b) \equiv \mod{p^k}
\begin{cases}
(n+a-1)^{p^k} \cdot T(K_{n -p^k};a,b) & n > p^k\\
(a-1)^{p^k-1} & n = p^k
\end{cases}
$$
\end{enumerate}
\end{conjecture}

\subsection{Outline of the paper}

In Section \ref{se:sb} we present the Specker-Blatter Theorem.
In Section \ref{se:matching} we prove Theorem \ref{th:3},
and in
In Section \ref{se:tutte} we prove Theorem \ref{th:1} and \ref{th:4}.
In Section \ref{se:msol} we prove Theorem \ref{th:2}.
Finally, in Section \ref{se:conclu} we present our conclusions and suggestions
for further research.

\section{The Specker-Blatter Theorem}
\label{se:sb}

The Specker-Blatter Theorem from 1980 is the first
application of Logic to Combinatorial Counting. At the time the theorem was hardly noticed,
mostly due to its unlucky placement for publication, \cite{pr:BlatterSpecker81,pr:BlatterSpecker84,ar:Specker88}.
An easily accessible place to find proofs and a survey of further developments is
\cite{ar:FischerKotekMakowsky11}.\vspace*{-2mm}

\subsection{Counting graphs: The density function}

A  {\em labeled graph $G$} with $n$ vertices will have $V(G) =\{0,1, \ldots , n-1\} =[n]$.
There are $gr(n)=2^{n \choose 2} = 2^{\frac{n(n-1)}{2}}$
many graphs with $n$ vertices.

A graph property is a class of finite graphs closed under graph isomorphisms.
For a graph property $\mathcal{P}$ denote by $\mathcal{P}^n$
the graphs with $n$ vertices in
$\mathcal{P}$, and by
$d_{\mathcal{P}}(n) = \mid \mathcal{P}^n \mid$,
the number of graphs $G$ with $V(G)=[n]$ which are in $\mathcal{P}$.
$d_{\mathcal{P}}(n)$ is called the {\em  density function} of $\mathcal{P}$.
The density function counts {\em labeled} graphs.
Let $G$ consist of the vertices $[n]$ and one single edge which is not a loop.
There is,$\,$up to isomorphisms (disregarding the labels) one such graph,$\,$but there are $n(n-1)/2$ such labeled$\,$graphs.

A graph property $\mathcal{P}$ is hereditary if it is closed under induced subgraphs.
$\mathcal{P}$ is monotone if it is closed under (not necessarily induced) subgraphs.
If
$\mathcal{P}$ is hereditary or monotone, the density function $d_{\mathcal{P}}(n)$  of
$\mathcal{P}$ is also called the
{\em  speed}  of
$\mathcal{P}$, since it is an {\em ultimately monotone increasing function},
cf. \cite{ar:BaloghBollobasWeinreich2000,ar:BaloghBollobasWeinreich2001,balogh2002measures}.
Studying the possible growth rate of the speed of a graph property was initiated in 1994
by E. Scheinerman and J. Zito in \cite{scheinerman1994size}.
E. Specker and C. Blatter already in 1980 studied under what conditions the density function of graph properties
satisfy recurrence relations.
The definition of the density function can be extended to relational structures of
any vocabulary $\tau$. However, we will consider only vocabularies $\tau$ without function symbols.
\begin{examples}
\begin{enumerate}[(i)]
\itemsep=0.9pt
\item
If $\mathcal{P} = Graphs$  consists of all  simple graphs,
$$
d_{Graphs}(n) = 2^{n \choose 2}
$$
In the unlabeled case the function is rather complicated.
\item
If $\mathcal{P} = LinOrd$ consists of all linear orders.\vspace*{-1mm}
$$
d_{LinOrd}(n) =  n!
$$
In the unlabeled case we have  the constant function
with value $1$.
\item
If $\mathcal{P} = SqGrids$ consists of all square grids,\vspace*{-1mm}
$$
d_{SqGrids}(n) =
\begin{cases}
\frac{n!}{4} & \mbox{ if } n = m^2 \\
0 & \mbox{ else }
\end{cases}
$$
In the unlabeled case we have $1$ instead of $n!$ in the above expression.
\end{enumerate}
\end{examples}

\vfil\eject

For some graph properties $\mathcal{P}$
the density functions satisfies a linear recurrence relation over $\Z$.
However, this is not always the case.

\begin{lemma}[Folklore]
\label{le:rec-growth}
Let $f: \mathbb{Z} \rightarrow \mathbb{Z}$ a function which satisfies
a linear recurrence relation
$$f(n+1) = \sum_{i=0}^k a_i f(n-i)$$
over $\mathbb{Z}$.
Then
there is a constant $c \in \mathbb{Z}$ such that $f(n) \leq 2^{cn}$.
\end{lemma}

\begin{corollary}
For $\mathcal{C} \in \{ Graphs, LinOrd, SqGrids\}$,
$d_{\mathcal{C}}(n)$ does not satisfy a linear recurrence over
$\mathbb{Z}$.
\end{corollary}

\subsection{Modular counting}

Let $\mu \in \N$.
\begin{obs}
For every $\mu \in \mathbb{N}$ and for large enough $n$ we have
$n! =0 \pmod{\mu}$
\end{obs}
Hence, for $n \geq N(\mu)$ we have
$$
d_{LinOrd}(n+1) = d_{LinOrd}(n) \pmod{\mu}\vspace*{-2mm}
$$
and\vspace*{-1mm}
$$
d_{SqGrid}(n+1) = d_{SqGrid}(n) \pmod{\mu}\vspace*{1mm}
$$
We say that a function $f(n)$ satisfies a {\em trivial modular recurrence}
if for every $\mu$ there exists $N_{\mu}$ such that
if $n>N_{\mu}$ then $f(n)\equiv 0 \pmod{\mu}$.
Clearly,
the two examples above satisfy
trivial modular recurrences.

\begin{obs}
$f(n)$ satisfies a trivial modular recurrence
iff
there exist functions $g(n), h(n)$ with $g(n)$ tending to infinity
such that $f(n) = g(n)! \cdot h(n)$.
\end{obs}
In other words,
trivial modular recurrences are always caused by some factor which is a factorial.

\medskip
It is sometimes more intuitive to say of integer sequences with values in $[m]$ that they
are ultimately periodic rather than to talk about modular linear recurrence.

\begin{proposition}[Folklore]
Let $a_n$ be an integer sequence.
\\
$a_n$ satisfies a linear recurrence relation modulo $\mu$
iff
$a_n$ is ultimately periodic modulo $\mu$.
\end{proposition}

The following are two instructive examples.
First, we look at the class of all graphs.
\begin{example}
\label{ex:1}
The density function for all graphs is given by
$$
d_{Graphs}(n+1) =
2^{n+1 \choose 2} =
2^{n \choose 2} \cdot 2^n.
$$
Therefore\vspace*{-2mm}
$$
d_{Graphs}(n+m) =
d_{Graphs}(n)  \cdot \prod_{i=0}^{m-1} 2^{n+i} =
d_{Graphs}(n)  \cdot 2^{nm} \cdot \prod_{i=0}^{m-1} 2^{i}
$$
As $a^{p-1}  =a \pmod{p}$ (Fermat's Little Theorem) we get
with $a=2^n$ and $m=p$ a prime
$$
d_{Graphs}(n+p) =
d_{Graphs}(n)  \cdot \prod_{i=0}^{p-1} 2^{i} \pmod{p}
$$
This is a non-trivial recurrence for $\mu=p$ a prime.
It is also different for distinct primes $p$ and $p'$,
In other words, the existence of the modular recurrence is non-uniform in $p$.
\end{example}

The second example is the class of graphs
$EQ_2CLIQUE$
which consists of the graphs which are the disjoint unions of two equal-sized cliques.

\begin{example}
\label{ex:2}
For density function $d_{EQ_2CLIQUE}(n)$ we have
$$
d_{EQ_2CLIQUE}(n) = b_2(n)=
\begin{cases}
\frac{1}{2} {2m \choose m} & \mbox{ for } n=2m
\\
0 & \mbox{ else }
\end{cases}
$$
The factor $\frac{1}{2}$ is there because we cannot distinguish
the choice of the first clique from the choice of its complement.
\end{example}

The function $b_2(n)$ was studied by
F. E. A. Lucas (1842-1891)
in 1878, but not published at the time.
It was found in his notes in the National Archive of France.
A proof may be reconstructed from the hints in \cite[Exercise 5.61]{bk:GrahamKnuthPatashnik94}.
\begin{proposition}[Lucas, 1878]
For every $n$ which is not a power of $2$, we have
$b_2(n)\equiv 0 \pmod{2}$, and for every $n$ which is a power
of $2$ we have $b_2(n)\equiv 1 \pmod{2}$.
\\
In particular, $b_2(n)$ is not ultimately periodic modulo $2$.
\end{proposition}
We conclude that
$d_{EQ_2CLIQUE}(n)$ is not ultimately periodic modulo $2$.

Finally, here is an example, were the precise counting is known. It is originally due to
Redfield, \cite{ar:Redfield27} and was rediscovered by R.C. Read and G. Polya, \cite{ar:Read59,ar:Read60,bk:PolyaRead}.
\begin{example}
\label{ex:3}
Let $\mathcal{R}_d$ be the class of regular graphs of degree $d$
and $d_{{\mathcal R}_d}(n)$ its density function.
\end{example}
It is not at all obvious that $d_{{\mathcal R}_d}(n)$ is ultimately periodic modulo $\mu$, but it follows
from the Specker-Blatter Theorem \ref{th:sb} below that, indeed, it is.
\begin{proposition}[J.H. Redfield, 1927]
For $d=3$ we have
$d_{{\mathcal R}_3}(2n+1) = 0$ and
$$
d_{{\mathcal R}_3}(2n) =\frac{(2n)!}{6^n}?\sum_{j,k}
\frac{(-1)^j(6k-2j)!6^j}{(3k-j)!(2k-j)!(n-k)!}48^k\sum_i\frac{(-1)^i j!}{(j-2i)! i!}
$$
\end{proposition}
An accessible proof can be found in
\cite[page 187]{bk:HararyPalmer}.

\subsection{$\MSOL$-definable graph properties}

We now look at the density functions of graph properties $\mathcal{P}$
which are definable in $\MSOL$.

\begin{theorem}[Specker-Blatter Theorem]
\label{th:sb}
Let $\mathcal{P}$
be a graph property which is $MSOL$-definable
and let
$d_{\mathcal{P}}(n)$ be its density function.
\begin{itemize}
\itemsep=0.95pt
\item
$d_{\mathcal{P}}(n)$ satisfies modular recurrence relations for
each $\mu \in \N$, hence it is ultimately periodic modulo $\mu$.
\item
This remains true for vocabularies $\tau$ with several binary edge relations and unary predicates
on the vertices.
\end{itemize}
\end{theorem}

The Specker-Blatter Theorem does not hold
if one allows quaternary relations in $\tau$,
\cite{ar:Fischer02,ar:Specker2005}.
\begin{theorem}[E. Fischer, 2003]
Let $\tau_0$ consist of one quaternary relation.
There is a class of
$FOL(\tau_0)$-definable  $\tau_0$-structures $\mathcal{F}$
such that
$d_{\mathcal{F}}$ is not ultimately periodic modulo $2$.
\end{theorem}
The proof consists of a very clever encoding of
$EQ_2CLIQUE$ using the quaternary relation.

The restriction to binary relations is not needed if the graph property $\mathcal{P}$ contains only graphs of bounded degree,
see \cite{pr:FischerMakowsky03}. The same works for classes of $\tau$-structures of bounded degree, where
the degree is defined via the Gaifman graph of the structures.

\medskip
The Gaifman graph of a $\tau$-structure $\aA$ is the (undirected, loop-free) graph $G_{\aA}$ with
vertex set $A$, the universe of $\aA$ and an edge between two distinct vertices $a, b in A$ iff there exists an $R  \in \tau$
and a tuple $(a_1, . . , a_r) \in R^A$ such that $a, b \in \{a_1, . . , a_r\}$, cf. \cite{bk:Libkin2004}.

\section{The matching polynomial of the complete graph}
\label{se:matching}

Here we prove
\begin{quote}
{\bf Theorem \ref{th:3}}:
For all $a,b, \mu \in \N$.
$\bar{M}(K_n;a,b)$
is ultimately periodic modulo $\mu$.
\\
In particular $M(K_n;a)$
is ultimately periodic modulo $\mu$.
\end{quote}
\begin{proof}
First we prove it for $M(K_n;a)$ and note first that $M(K_n;a)$
is of the form
$$
M(K_n;a) =
P_{\psi}(K_n,a) = \sum_{F \subseteq E(G): \psi(E,F)} a^{|F|}
$$
where $\psi(E,F)$ says that $E$ is the edge relation of $K_n$ and $F \subseteq E$ is a set of independent edges of $K_n$.
We interpret $a^{F}$ as the set of functions $f: F \rightarrow [a]$.
Each $f$ induces a  partition of $F$ with
$$F_f(i) = \{ e \in F: f(e) =i \}.$$
We have
$$
a^{|F|} = |\{ f: F \rightarrow [a] \}|.
$$
Next
we look at the density function of
$$
g_{a}(n) = |\{ U_1, \ldots U_a \subseteq [n]^2, F \subseteq [n]^2 : \phi_1(\bar{U},F), \psi(E,F) \} |
$$
where
$\phi_1(\bar{U})$ says: $U_1, \ldots, U_{a}$ partition $F$ and
$\psi(E,F)$ says that $E$ is the edge relation of $K_n$ and $F \subseteq E$ is a set of independent edges of $K_n$.
$g_{a}(n)$ encodes the computation of $M(K_n;a)$.
\begin{claim}
$$
g_{a}(n) = P_{\psi}(K_n,a) = \sum_{F \subseteq E(G): \psi(E,F)} a^{|F|} =M(K_n;a)
$$
\end{claim}
All the relation symbols of $\phi_1$ and $\psi$
are binary, therefore we can apply the Specker-Blatter Theorem and conclude that
$g_{a}(n)$ is ultimately periodic modulo every $\mu \in \N$, and so is $M(K_n;a)$.

\medskip
The proof for $\bar{M}(K_n;a,b)$ is similar.
$$
\bar{M}(K_n;a,b) =
\sum_{k=0}^{\lfloor n/2 \rfloor} (X)^k m_k(G) Y^{n-2k} =
\sum_{F \subseteq E(G): \psi(E,F)}
(a)^{|F|} m_k(G) b^{|dF|}
$$
where $dF$ is the set of vertices in $[n]$ not covered by $F$.
If $F$ has $k$ edges, then $dF$ has $n-2k$ vertices.
\end{proof}

\vspace*{-3mm}
\section{The Tutte polynomial of a complete graph}
\label{se:tutte}

Now we prove
\begin{quote}
{\bf Theorem \ref{th:1}}
\\
For every $a,b, \mu \in\N^{+}$ with $a,b>1$,
$\gcd(a-1, \mu )=1$ and $\gcd(b-1, \mu )=1$,
$T_n(a,b)=T(K_{n},a,b)$ is ultimately periodic modulo $\mu$.
\end{quote}

\begin{proof}
We first rewrite the Tutte polynomial as
$$
T(G;X,Y) =
\frac{1}{(Y-1)^{|V(G)|} \cdot (X-1)^{\kappa(G))}}
\sum_{A \subseteq E(G)} (X-1)^{\kappa(A)} \cdot (Y-1)^{\kappa(A) + |A|}
$$
and put $G =K_n$. Now $|V(K_n)|=n$ and $|E(K_n)|={n \choose  2}$, and  $\kappa(K_n)=1$.
$$
T(K_n;X,Y) =
\frac{1}{(Y-1)^{n} \cdot (X-1)}
\sum_{A \subseteq E(G)} (X-1)^{\kappa(A)} \cdot (Y-1)^{\kappa(A) + |A|}.
$$
\eject
Like in the proof Theorem \ref{th:3}
we interpret
$a^{A}$ as the set of functions $f: A \rightarrow [a]$.
Each $f$ induces a  partition of $A$ with
$$A_f(i) = \{ e \in A: f(e) =i \}.$$
We have
$$
a^{|A|} = |\{ f: A \rightarrow [a] \}|.
$$
In the case of
$a^{\kappa(A)}$ the sets $A_f(i)$ have to be $A$-closed, in order to partition the connected components
of the spanning induced by the set of edges of $A$.
Next look at the function
$$
f_{a,b}(n)= (b-1)^{n} \cdot (a-1) \cdot T(K_n, a,b)
=
\sum_{A \subseteq E(G)} (a-1)^{\kappa(A)} \cdot (b-1)^{\kappa(A)} \cdot (b-1)^{|A|}
$$
as a  $\MSOL$-definable density function:
\begin{gather}
f_{a,b}(n) = |\{\bar{U}, \bar{R}, \bar{S} \subseteq [n]:
\phi_1(\bar{U},A), \phi_2(\bar{R}, A), \phi_3(\bar{S}), A \subseteq E \} |
\notag
\end{gather}
where for
$\bar{U} = (U_1, \ldots, U_{a-1})$ and
$\bar{R} = (R_1, \ldots, R_{b-1})$ are unary relations and
$\bar{S} = (S_1, \ldots, S_{b-1})$ are binary, and
\begin{enumerate}[(i)]
\item each $U_1, \ldots, U_{a-1}, R_1, \ldots, R_{b-1}, S_1, \ldots, S_{b-1} \subseteq [n]$;
\item  $\phi_1(\bar{U}, A)$ says: $U_1, \ldots, U_{a-1}$
partitions the connected components of $G[A]$,
and each $U_i$ is a disjoint union of connected components of the graph $G= ([n], A)$;
\item
$\phi_2(\bar{R}, A)$ says: $R_1, \ldots, R_{b-1}$ also
partitions the connected components of $G = ( [n], A)$,
\item and
$\phi_3(\bar{S})$ says: $S_1, \ldots, S_{b-1}$ partitions $A \subseteq [n]^2$.
\end{enumerate}
All the formulas
$\phi_1(\bar{U}, A), \phi_2(\bar{R}, A), \phi_3(\bar{S})$
are in $\MSOL$
and
contain only unary and binary relation symbols.
It follows by Theorem \ref{th:sb} that $f_{a,b}(n)$ is
ultimately periodic modulo every $\mu \in \N$.
\begin{claim}
$$
f_{a,b}(n) = \sum_{A \subseteq E(G)} (a-1)^{\kappa(A))} \cdot (b-1)^{\kappa(A)} \cdot (b-1)^{|A|}
$$
\end{claim}

Now we need a lemma.
\begin{lemma}
\label{lemma}
Let $d_1(n), d_2(n)$ be integer functions.
\begin{enumerate}[(i)]
\item
Let $c , \mu \in \N^+ $.
\\
Assume $c \cdot d_1(n)$ is ultimately periodic modulo $\mu$ and $\gcd(c, \mu)=1$.
\\
Then $d_1(n)$ is ultimately periodic modulo $\mu$.
\item
Let $t, \mu \in \N^+$ with $\gcd(t,\mu)=1$.
\\
Assume that $t^n \cdot d_2(n)$ is ultimately periodic modulo $\mu$.
\\
Then $d_2(n)$ is ultimately periodic modulo $\mu$.
\end{enumerate}
\end{lemma}
\begin{proof}
(i) is left to the reader.
\\
(ii):
Since $t$ and $\mu$ are relatively prime, $t$ has a multiplicative
inverse $g$ modulo $\mu$: $g \cdot t \equiv 1 (\mod{\mu})$.
The product of two ultimately periodic sequence is ultimately
periodic, hence this is true for the product of $g$ and  $t^n \cdot d_2(n)$.
But we have
$$g^n \cdot t^n \cdot d_2(n) \equiv d_2(n)(\mod{\mu}),$$ hence $d_2(n)$
is ultimately periodic.
\end{proof}
To complete the proof of Theorem \ref{th:1} we note that $\kappa(E) =1$  and $|V|=n$,
and we put:
$$
d_1(n) = T(K_n, a,b) \mbox{  and  } d_2(n) = (b-1) \cdot T(K_n, a,b)
$$
We have $f_{a,b}(n)= (a-1)^n \cdot d_2(n)$ is ultimately periodic modulo $\mu$
for $\gcd((a-1),\mu)=1$.
By the Lemma  \ref{lemma}(ii) we have that $d_2(n)$ ultimately periodic modulo $\mu$
for $\gcd((a-1),\mu)=1$.
By the Lemma  \ref{lemma}(i) we have that $d_1(n)$ ultimately periodic modulo $\mu$
for $\gcd((b-1),\mu)=1$.
\end{proof}

The proofs of Theorem \ref{th:4}  for $S(G;a,b,c)$ and $C(G;a,b,c)$ are similar as for  $f_{a,b}(n)$.
\\
For $\xi(G;a,b,c)$ we note that $\xi((G;X,Y,Z) = C(G;X,Y, \frac{Z}{XY} +1)$ from Proposition \ref{trinks}.
\\
To be able to use now $C(G;a,b, \frac{c}{ab} +1)$ one requires that $ab$ divides $c$ so that
$\frac{c}{ab}$ is a non-negative integer.

\section{Theorem \ref{th:2} and its limitations}
\label{se:msol}

The logic $\CMSOL$ is the extension of $\MSOL$ by adding {\em modular counting quantifiers}.
Let $\phi(x)$ be a formula with free variable $x$.
A modular counting quantifier $C_{\mu, k}x \phi(x)$
says that there are, modulo $\mu$, $k$-many elements satisfying $\phi(x)$.
Using Ehrenfeucht-Fra\"iss\'e games for $\MSOL$ one can show that
$C_{\mu, k}x \phi(x)$ is not expressible in $\MSOL$. However,
for $\phi(x)$ a formula in $\SOL$ the formula
$C_{\mu, k}x \phi(x)$ it is expressible in $\SOL$.
We say that there is an equivalence relation $E$
on the set defined by $\phi(x)$ which has exactly one equivalence class of size $k$ and all the other
non-empty equivalence classes are of size $\mu$. There are other ways of defining $C_{\mu, k}x \phi(x)$
in $\SOL$. Instead of the existence of one binary relation we can also assert that there are $\mu$ disjoint subsets
of equal size
of the set defined by $\phi(x)$
and the complement of their union has size $k$.
But then we need the existence of a binary relation which expresses
that the unary predicates are of equal size.

\medskip
$\CMSOL$ is like $\MSOL$ but modular counting quantifiers are allowed.
The syntax of $\CMSOL$ is obtained from the syntax of $\MSOL$ by allowing
also quantification with $C_{\mu, k}x$. The meaning function of $\MSOL$ is then
naturally extended to $\CMSOL$.
$\CGMSOL$ is defined analogously by adding modular counting quantifiers to $\GMSOL$.

It was shown in \cite{pr:FischerMakowsky03,ar:FischerKotekMakowsky11}
that
the Specker-Blatter Theorem also holds for $\CMSOL$ for relational structures with
relations of arity at most two, and for $\CGMSOL$ graphs.

\subsection{$\CMSOL$-definable graph polynomials}

We now discuss how Theorems \ref{th:1} and \ref{th:3} can be extended.

\medskip
A univariate graph polynomial of the form
$$
P_{\phi}(G;X) =\sum_{A \subseteq V(G): \phi(E,A)} X^{|A|}
$$
is $\CMSOL$ definable if $\phi(A,E)$
is a $\CMSOL$-formula in the language of graphs with an additional predicate for $A$.

\medskip
A univariate graph polynomial of the form
$$
P_{\psi}(G;X) =\sum_{F \subseteq E(G): \psi(E,F)} X^{|F|}
$$
is $\CGMSOL$ definable if
$\psi(F,E)$
is a $\CGMSOL$-formula in the language of graphs with an additional predicate for $A$
or for $F \subseteq E$.
The independence polynomial and the clique polynomial
are of the first form.
The matching polynomials $\alpha$ and $M$ are of the second form.

\medskip
Analyzing the proof of Theorem \ref{th:3} immediate gives Theorem \ref{th:2}:
\begin{quote}
If $P$ is a graph polynomial definable in $\CGMSOL_{graph}$ without order,
then for every $\mu \in \N$ and every $a,b \in \Z$\\
the sequence $P(K_n, a,b)$ is ultimately periodic modulo $\mu$.
\end{quote}
The proof of Theorem \ref{th:1} also works for the graph polynomials listed in Theorem \ref{th:4}.
Unfortunately, we have not found interesting applications of Theorem \ref{th:2}.
In many cases $P(K_n, a,b)$ can be shown directly to
satisfy linear recurrence relations over $\Z$.
\begin{examples}
\begin{enumerate}[(i)]
\itemsep=0.9pt
\item
If $\mathcal{P}$ is a graph property which does not contain any complete graph and is
definable in $\CMSOL$ for graphs
and $\phi_{\cP}(A,E)$ says that $A$ induces a graph $G[A] \in \cP$
then
$$
P_{\phi_{\cP}}(G;X) =\sum_{A \subseteq V(G): \phi_{\cP}(E,A)} X^{|A|}
$$
trivializes for $P_{\phi_{\cP}}(K_n;X) $.
\item
If instead, $\phi^{'}_{\cP}(A,E')$ says that $(A, E')$ is a subgraph  $\cP$
$P_{\phi^{'}_{\cP}}(K_n;X)$ may be non-trivial.
\item
The domination polynomial $D(G;X)$ is obtained by taking
$\phi_{dom}(A,E)$ which says that $A$ is a dominating set for $E$.
In \cite{alikhani2009introduction,oboudi2016roots}
it is shown that $D(K_{n+1};X) = D(K_n;X)(X+1) +X$.
For $X =a$ this is a linear recurrence relation, hence Theorem \ref{th:2}
gives nothing new.
\item
The univariate interlace polynomial $q(G;X)$ from \cite{arratia2004interlace} is $\CMSOL$-definable, as shown in
\cite{courcelle2008multivariate}. But  $q(K_n;X) = 2^{n-1}X$, hence $q(K_{n+1};X) = 2X q(K_n;X)$.
Again Theorem \ref{th:2} gives nothing new.
\end{enumerate}
\end{examples}

\begin{problem}
Find  more graph polynomials $P(G; \bar{a})$  where for $\bar{a} \in \N^k, \mu \in \N$ the
sequence $P(K_n; \bar{a})$ is not obviously ultimately periodic modulo $\mu$.
\end{problem}

\begin{problem}
Find interesting cases of graph polynomials  where Theorem \ref{th:2} gives a non-trivial result.
\end{problem}

We have seen that Theorem \ref{th:2} also holds for multivariate polynomials.
The Tutte polynomial is not of this form, because of the term $(X-1)^{\kappa(S)}$.
It is of this form in the language of ordered graphs.
However, the Specker-Blatter Theorem formulated for ordered structures trivializes, because
the ordering adds a factor of the form $n!$, hence the sequence satisfies a trivial modular recurrence.

\subsection{Why complete graphs?}

Complete  labeled graphs on $n$ vertices are definable in $\FOL$ in the empty vocabulary
and are unique, not only up to isomorphisms.
The same is true about the empty graph.

\begin{proposition}
Assume $\phi(\bar{x},\bar{y})$ is a formula of $\CMSOL$ over the empty vocabulary which
defines a unique edge relation  $E_{\phi}$ on $k$-tuples of a set $[n]$.
Then the graph $G = ([n]^k,  E_{\phi})$ is either a complete graph or the empty graph.
\end{proposition}
The way we used the Specker-Blatter Theorem  in the proof of Theorem \ref{th:2} did require that the function
$$
g_{a}(n) = |\{ U_1, \ldots U_a \subseteq [n], F \subseteq [n]^2 : \phi_1(\bar{U},F), \psi(F) \} |
$$
evaluates the matching polynomial at $a$ for $G = K_n$.
If instead of $G = K_n$ we use some $G = ([n]^k,  E_{\phi})$ with $k=1$ we can modify $g_{a}(n)$
and still apply the Specker-Blatter Theorem to it, but the modified version $g_{a, \phi}(n)$ will
not evaluate the matching polynomial anymore.
Let $sp(\phi, n)$ be the number of labeled graphs of the form
$G = ([n],  E_{\phi})$.
Even if all the graphs
$G_n = ([n],  E_{\phi})$ are isomorphic, $g_{a, \phi}(n)$
computes
$$
g_{a, \phi}(n)=
M(G_n;a) \cdot sp(\phi, n)
$$
which will be ultimately periodic modulo $\mu$. However, this does not suffice to conclude
that
$M(G_n;a)$ is ultimately periodic modulo $\mu$.

\medskip
{\bf Facit:} Our proofs of Theorems
\ref{th:1},
\ref{th:2} and
\ref{th:3}
only work for $G=K_n$ or $G=\bar{K_n}$.

\section{Conclusions}
\label{se:conclu}

Inspired by the paper of A.P. Mani and R.J. Stones \cite{ar:ManiStones2016}
we have examined modular recurrences for the matching and the Tutte polynomial
of a complete graph. We have noted that the existence of modular recurrence relation
follows from the Specker-Blatter Theorem, without explicitly describing the
exact modular recurrence. The conjectures of A.P. Mani and R.J. Stones
are more ambitious, as they give a precise statement about how these modular recurrences
look in the case of the Tutte polynomial.
We also noted that our approach via the Specker-Blatter Theorem \ref{th:sb}
works for other graph polynomials definable in variants of $\MSOL$.
Most strikingly, it works for the trivariate edge elimination polynomial $\xi(G;X,Y,Z)$,
which is the most general edge elimination polynomial, and generalizes both the matching polynomials
and the Tutte polynomial.

Our result suggests that the use of the Specker-Batter Theorem to
$\CMSOL$-definable ($\CGMSOL$-definable) graph polynomials should be further
investigated.
In \cite{fischer2003specker,ar:FischerKotekMakowsky11} a logic-free version of the Specker-Batter Theorem is
discussed, where instead of the graph property $\cP$  being definable in $\CMSOL$
one only requires that $\cP$ has finite Specker rank.
Here we note that there uncountably many graph properties of finite Specker rank,
but there only countably many $\CMSOL$-definable graph properties.
$\CMSOL$-definable graph properties have always finite Specker rank, but the upper bound on the Specker rank
which follows from $\CMSOL$-definability is very often exponentially bigger than the true Specker rank.
On the other hand the parameters which characterize the ultimate periodicity depend only polynomially
on the Specker rank, hence
the structure of the modular recurrences can be more precisely described
using the Specker rank.

\subsection*{Acknowledgments}

The material of this paper was first presented at the
Colloquium Logicum, September 10-12, 2016 held in Hamburg, Germany.
An updated version was given at the Colloquium in honor of B. Trakhtenbrot's centenary,
held in Tel Aviv (and via ZOOM) in October 2021.
The second author also presented a survey of the Specker-Blatter Theorem and its
various extensions at the Colloquium held at ETHZ in honor of E. Specker's centenary in February 2020
in Zurich, Switzerland.
Both E. Specker and C. Blatter were my (the second author's) teachers in 1967-1973. They both shaped my
mathematical and logical views in a lasting way. Boris Abramovich did the same after we
both started to live in Israel after 1980. He gave his unrelenting moral support for my own work.
Thanks also to the unknown referee for his careful reading and his resulting suggestions, which helped improving the
paper. Thanks also to Peter Tittmann for useful comments.


\section{\label{sec:level1}First-level heading}
Many stochastic processes such as radioactive decay or the inter-arrival times of new bitcoin blocks are memoryless, i.e. the probability of an event to occur is constant in time.  Fluid flows on the other hand are governed by deterministic equations, and the flow state at an instance in time in principle fully determines the entire future evolution. It may therefore come as a surprise that the transition to turbulence in a large number of flows is governed by memoryless processes. This seemingly contradictory state of affairs is reconciled by the sensitive dependence on initial conditions --- a defining property of deterministic chaos that equally applies to turbulent fluid motion. The associated exponential amplification of even minute differences makes long term forecasts of the dynamics impossible in practice, and effectively erases the flow's memory.

It is precisely this memoryless property that has paved the way for a possible solution to the century old puzzle surrounding the nature of the transition to turbulence.  More specifically, the transition in Couette, pipe and related shear flows is characterized by the co-existence of laminar and turbulent regions. Such patches of turbulence, although frequently long lived, are transient, and their decay is abrupt and unpredictable --- more precisely, memoryless. The underlying dynamical-systems model is that of a chaotic saddle in phase space \cite{eckhardt2007turbulence}, and the laminar state  corresponds to an attracting fixed point, where the escape rate from the former is constant in time \cite{tel2008chaotic}. This is probably the best established property of transitional turbulence in subcritical shear flows, and a large body of literature (e.g.\citep{bottin1998discontinuous,faisst2004sensitive,moehlis2004low,peixinho2006decay,hof2006finite,willis2007critical,hof2008repeller,schneider2008lifetime,avila2010transient,borrero2010transient,linkmann2015sudden,shimizu2019exponential}) has been dedicated to the dynamical-systems aspects, including studies of the underlying exact coherent structures\cite{hof2004experimental,kerswell2005recent,eckhardt2007turbulence,kawahara2012significance}, the boundary crisis \cite{van2011homoclinic,kreilos2012periodic,avila2013streamwise,ritter2016emergence,budanur2019geometry} that gives rise to transient chaos, and the resulting edge state that mediates the eventual escape from turbulence \cite{itano2001dynamics,schneider2007turbulence}.
Corresponding transient dynamics have equally been observed for turbulence in astrophysical accretion discs\cite{Rempel2010}, in nonlinear fibre optics\cite{Falkovich2013}, for the dynamo transition in magnetohydrodynamic turbulence\cite{Rempel2009}, turbulence in pulsatile\cite{Xu2017} (e.g. cardiovascular) flows and surprisingly even extend to the generic case of isotropic turbulence\cite{linkmann2015sudden}.

In addition to the temporal dynamics, also the spatial proliferation, i.e. the nucleation of new patches via a contact process (called `puff splitting' \cite{moxey2010distinct} in pipes), turned out to be memoryless \cite{avila2011onset}. These combined insights have not only allowed to determine the critical point for the onset of sustained turbulence \cite{avila2011onset}, but also led to a direct analogy to non-equilibrium phase transitions and in particular to directed percolation (DP): individual patches of turbulence can be interpreted as active sites which just like in models of DP can either decay or infect their neighbours, and the probability for either event to take place does not depend on time. Based on these insights, DP is believed to be the standard route to turbulence for this entire class of flows \cite{shi2013scale,lemoult2016directed,sano2016universal,chantry2017universal,hiruta2020subcritical,peinkePRX}. Moreover these concepts have been applied to the dynamics of active biological matter \cite{Doostmohammadi2017} as well as to the Leidenfrost effect\cite{Chantelot2021} in soft matter physics where starting from observations of Markovian structures (analogous to puffs) DP exponents have been measured.

At first sight, channel flow perfectly complies with all of the above criteria. In the transitional regime, turbulence appears in the form of stripes, and simulations in periodic boxes confirmed the memoryless nature of decay and splitting of individual stripes \cite{shimizu2019exponential,gome2020statistical,tuckerman2020patterns}. Like in other shear flows \cite{avila2011onset,shi2013scale}, the critical point has been approximated by comparing mean life and splitting times \cite{gome2020statistical}.
Additionally, an experimental study \cite{sano2016universal} suggested that the onset of sustained turbulence falls into the DP universality class seemingly explaining the nature of this transition. A closer look at the latter study however reveals that the time scales resolved ($\sim 2400$ advective time units) were more than 3 orders of magnitude shorter than those proposed in the DP analogy for pipe and Couette flow. Moreover the spot like flow structures formed during the short time scales of this experiment did not correspond to the large scale stripe patterns that are characteristic for transitional channel flows (e.g. see \cite{xiong2015turbulent,shimizu2019bifurcations} and Fig. 1). Finally as shown by a number of recent studies \cite{tao2018extended,xiong2015turbulent,shimizu2019bifurcations,manneville2020transitional}, the actual critical point where turbulent stripes first become sustained in channel flow is substantially lower than the critical point proposed by \cite{sano2016universal} and consequently these latter experiments were carried out far from the onset of sustained turbulence.

A model study, attempted to remedy the shortcomings of these channel experiments by suggesting that at the actual (i.e. much lower) critical point the transition may fall into the 1+1D DP universality class and crosses over to 2+1 D DP further from critical   \cite{shimizu2019bifurcations,manneville2020transitional}.
On the other hand, a numerical study of nonlocalized stripes in periodic domains \cite{gome2020statistical} reported that the characteristic time scales for the decay and the nucleation of new stripes exceed $10^6$ advective time units at the critical point. While it is not clear if the decay and proliferation of fully localised stripes follow the same rules, such excessive time scales would put a characterisation of the nature of the transition practically beyond reach.

We will show in the following that
stripe localization a process intrinsic to transitional channel flow, fundamentally alters the proliferation and decay processes of turbulence. Unlike for other shear flows, stripes in channel flow become sustained via a deterministic growth process. Even close to the critical point, and over the longest observation times realizable, turbulent stripes are not memoryless, the critical point cannot be approximated in the standard way \cite{avila2011onset}, the common edge tracking routine does not converge and the aforementioned mapping of turbulent stripes to DP sites does not hold.
\newline

\begin{figure}[h]
 \includegraphics[width=0.7\textwidth]{./figure_1}
 \caption{{\bf Transition from shrinking to expanding turbulent stripes.} The mean growth rate of turbulent regions is plotted as a function of Re. The first 4 points are fitted by a straight line to estimate the Re at which the growth rate first becomes positive, which is $Re \sim 650$. Above this Re, individual stripes grow, while below, they eventually decay, as seen in the insets showing a sequence of flow visualisation images. }
 \label{fig:figure1}
\end{figure}

Experiments are carried out in a large aspect ratio channel, considerably exceeding domain sizes of most earlier studies. Specifically, the channel's dimensions are $(Lx; Ly; Lz) = (4000h; 2h; 490h)$, where Lx, Ly and Lz are lengths in the streamwise, wall-normal and spanwise directions respectively. Quantities are non-dimensionalised with the length scale $h$ (half-gap) and 1.5 times the bulk velocity $u_{b}$  (corresponding to the centre-line velocity in case of laminar flow). In keeping with numerical and theoretical studies, we define the Reynolds number as $Re = 1.5 u_{b} h / \nu$, where $\nu$ is the kinematic viscosity. Unless perturbed, the flow stays laminar over the entire Reynolds number range investigated.

We used two protocols to obtain isolated stripes. Stripes were created either by perturbing the flow across a certain width ($\sim 20h$ in the spanwise direction) or a stripe was thus created at a higher flow-rate, followed by a quench to the target Re.  The latter mechanism was found to be more efficient at lower Re, at which the former creates stripes only sporadically. At $Re$ where both mechanisms work well, results are qualitatively the same regardless of the perturbation type, provided that the created stripes are of comparable size.

Turbulent stripes are investigated in the range $600 \leq Re \leq 900$. For each $Re$, around $1000$ individual stripes are generated near the channel entrance.

As they advect downstream, they are monitored by multiple cameras along the channel, until they eventually exit the channel. More details are available in the Supplementary Information. Note that in all the experiments, stripes remained at a minimum distance of $50 h$ from the lateral boundaries and stripe advection velocities and tip speeds remain unchanged up to this point indicating that stripes are unaffected by the side walls.
Overall more than $10,000$ stripes are tracked in order to obtain reliable statistics and the total observation time exceeds $5\times 10^7$ advective time units.  The resulting images are analysed to extract various quantities of interest. In particular, the area occupied by turbulence is determined as a function of time and used to compute the mean growth rate of turbulence, which is shown in  Fig.\ref{fig:figure1}. For $Re\lesssim650$, growth rates are negative, while above this Re, the growth rate becomes positive and continues to increase with $Re$. We propose that once individual stripes continuously grow also turbulence as a whole will eventually cease to decay and hence the critical point for channel flow can be estimated to be $Re_c\approx 650$.  In agreement with recent numerical studies  \citep{tao2018extended,xiong2015turbulent,shimizu2019bifurcations,manneville2020transitional} our experiments also confirm that the critical point is significantly lower than that reported by \cite{sano2016universal}.

\begin{figure}[h]
\centering
\begin{subfigure}{0.6\textwidth}
\includegraphics[width=\textwidth]{./figure_2a}
\caption{}
\label{fig:figure2a}
\end{subfigure}
\qquad
\begin{subfigure}{0.6\textwidth}
\includegraphics[width=\textwidth]{./figure_2b}
\caption{}
\label{fig:figure2b}
\end{subfigure}
\caption{{\bf Aging of stripes.} Panel (a) shows the survival probability of turbulent structures in pipe and channel flow. Non-localized stripes from DNS (blue dots), show the expected exponential drop in survival probabilities signifying a memoryless decay similar to puffs in pipe flow (grey symbols, data taken from \cite{avila2010transient}).  Natural, i.e. fully localized stripes in experiments (red circles and triangles) or DNS (blue squares) show aging,  i.e a decay process that is not memoryless, as seen by survival probabilities that are clearly not exponential. (b) Evolution of turbulent kinetic energy of stripes in direct numerical simulations. In moderate domains (grey curves) where stripes cannot localize the decays occur at unpredictable times.  Fully localised stripes (coloured curves) show a continuously decreasing energy and decays occur at comparable times. }
\label{fig:figure2}
\end{figure}

In order to characterize the decay of turbulence below the critical point, we next investigate stripe lifetimes at $Re=600$, i.e. $\sim 8\%$ below $Re_c$. Fig \ref{fig:figure2a} shows the probability that a stripe survives until a time t.
It is apparent from Fig. \ref{fig:figure2a} (red dots) that survival probabilities do not follow the exponential distribution that characterizes a memoryless process, and would equally be expected for the statistics of DP sites. For comparison we show puff lifetimes \cite{avila2010transient} for pipe flow at $Re=1900$, $\sim 7\%$ below $Re_c$ and at an even lower value of $Re=1860$, $\sim 9\%$ below critical, and both data sets (grey symbols) clearly follow exponential distributions characteristic for a chaotic saddle. In channel flow on the other hand, even upon a closer approach to critical ($Re=620$, $\sim 5\%$ below $Re_c$), stripe lifetimes (Fig. \ref{fig:figure2a}, red triangles) show no sign of memorylessness. In contrast to puffs, the decay rate of localized stripes  is not constant: No decays are found during a long initial phase, that grows as the critical point is approached ($\sim 900$ advective time units at $Re=600$, $\sim 1300$ advective time units at $Re=620$ ). Subsequently all stripes decay in a comparably short time interval (within $\sim 500$ advective time units). It is also apparent from Fig. \ref{fig:figure2a} that the tail of the respective distribution (red stars and red triangles) does not approach a constant slope but continues to steepen down to the lowest probabilities measurable. While initial transient periods are known for pipe flow\citep{hof2006finite}, these are however only of order $100$ advective time units and mark the period during which the puff forms ($t_0$). In our channel experiments the stripes are already fully formed from the start and subsequently undergo a slow aging process and as will be shown below, during this phase, the stripe recreation, or more precisely the streak production rate, slows down.

The shape of the lifetime distributions do not only fundamentally differ from those in pipe and Couette flows \citep{hof2006finite,shi2013scale}, but surprisingly also from computations of channel flow in small domains. In such computational boxes stripes cannot fully localize and are found to be memoryless \cite{shimizu2019exponential,gome2020statistical}. To verify these computational studies, we carried out direct numerical simulations in a long yet slender channel (30h; 2h; 100h) that is tilted by $45^{\circ}$, the same angle that stripes naturally assume with respect to the mean flow; see Supplementary Information for details. This method \cite{barkley2005computational} ensures that stripes, just like in previous numerical studies, connect via the periodic boundary, and hence localization is prohibited. Just like observed in earlier simulations, such periodic stripes have exponentially distributed lifetimes and are hence memoryless (blue
dots in Fig.~\ref{fig:figure2a}). However when the computational domain size is increased (400h; 2h; 400h) to allow stripes to fully localize and assume their natural shape, the memoryless nature disappears and stripes age, confirming the experimental observation. As shown in Fig. \ref{fig:figure2a} (blue squares), at $Re=620$,  localised stripes do not decay for the first $\sim 900$ advective time units, but all of them disappear within a subsequent interval of $\sim 500$ time units, in excellent qualitative agreement with our experiments.
In addition to whole stripes we can also partition stripes and determine the survival probability of stripe segments (i.e. a fixed number of streaks within a stripe), and in this case as well the corresponding distributions are not memoryless.

It should also be noted that for such aging processes, the time of decay and hence the lifetime depends on the energy threshold chosen. This becomes evident from the energy time-series of stripes, shown in Fig. \ref{fig:figure2b} for $Re=620$. Non-localised stripes (grey lines) show a well defined energy plateau for turbulence and decays can happen at any time with equal probability. In contrast, localised stripes (coloured lines) age, exhibiting a gradual energy decay. The time of the decay is not clearly defined and depends on the energy cut off chosen (as shown in Supplementary Fig.4) which however does not affect the qualitative shape of the distributions .
The lack of a clear energy plateau also makes standard bisection methods \citep{schneider2007turbulence,toh2003periodic}, commonly used to track the edge between laminar and turbulent flow, unsuitable.

In order to determine if closer to critical, stripe decays may be memoryless, we conducted experiments at $Re=640$. However here, virtually no stripe decays could be observed over the length of the channel. Even after increasing the number of investigated stripes to more than 1000, no decay was detected, despite the fact that overall stripes were shrinking in size. Hence, the aging time appears to have increased beyond the advection time through the channel. For memoryless decays, limitations in the domain length can simply be compensated by an increase in the number of cases investigated (e.g. as carried out by \cite{avila2011onset}), this is however not applicable if structures age like in the present case.

We next investigate the proliferation of stripes slightly above the critical point, at $Re=655$. As shown in Fig.~\ref{fig:figure3a}, stripes do not just passively advect downstream but at the same time move across the channel, i.e. in the spanwise direction. The spanwise motion is caused by the continuous creation of streaks at the downstream tip (see Fig.\ref{fig:figure3c} and movie in the Supplementary Information). Unlike for growth processes in other shear flows (e.g. puff nucleation in pipe flow \cite{avila2011onset} and in particular also unlike for stripe growth in Couette flow \cite{duguet2011stochastic}), here streaks are created at a constant rate (Fig.~\ref{fig:figure3d}) and this process is hence deterministic. This streak creation has recently been linked to a hydrodynamic instability particular to localized stripes in channel flow \cite{xiao2020growth}.
If the Reynolds number is decreased below critical, streaks are initially created at the same rate which results in an identical stripe migration speed (compare blue dashed lines in Fig. \ref{fig:figure3a} and \ref{fig:figure3b}). As long as this streak production mechanism is intact, stripes do not decay. However below critical, the production eventually slows down, and the tip speed decreases (the tip deviates from the blue dashed line in \ref{fig:figure3b}). Consequently the stripe shrinks and eventually dies.

\begin{figure}
\centering
\begin{subfigure}[h]{0.48\textwidth}
 \includegraphics[width=\textwidth]{./figure_3a}
 \caption{}
 \label{fig:figure3a}
\end{subfigure}
\hspace{1mm}
\begin{subfigure}[h]{0.48\textwidth}
 \includegraphics[width=\textwidth]{./figure_3b}
 \caption{}
 \label{fig:figure3b}
\end{subfigure}
\vspace{1mm}
\begin{subfigure}[h]{0.48\textwidth}
\includegraphics[width=\textwidth]{./figure_3c}
\caption{}
\label{fig:figure3c}
\end{subfigure}
\hspace{2mm}
\begin{subfigure}[h]{0.48\textwidth}
\includegraphics[width=\textwidth]{./figure_3d}
\caption{}
\label{fig:figure3d}
\end{subfigure}
 \caption{{\bf Deterministic stripe expansion and streak creation.} Panel (a) displays consecutive positions of a stripe at $Re=655$.  The downstream tip moves at a constant velocity diagonally across the channel as indicated by the dashed blue line, and laminar fluid is entrained a constant rate (b)  Below the critical point, at $Re=620$, the tip initially moves with the same velocity as the stripe at $Re=655$ but slows down in time and eventually the stripe decays. (c) Snapshots of the wall normal velocities in the region around the downstream tip of the stripe at $Re = 655$, obtained from DNS. showing a pair of streaks are created every $\sim 18$ advective time units. To more easily follow the addition of streaks, the streaks emerging in the first snapshot (t=0, arbitrary time origin) is shown with an arrow in all the subsequent snapshots (d) The time interval between the addition of two successive streaks shown as a function of time for $Re = 655$ }
\label{fig:figure3}
\end{figure}

\begin{figure}[h]
\centering
\begin{subfigure}[h]{0.5\textwidth}
\includegraphics[width=\textwidth]{./figure_4a}
\caption{}
\label{fig:figure4a}
\end{subfigure}
\qquad
\centering
\begin{subfigure}[h]{0.5\textwidth}
\includegraphics[width=\textwidth]{./figure_4b.eps}
\caption{}
\label{fig:figure4b}
\end{subfigure}
\caption{{\bf Age dependence of stripe splitting.} (a) Probability of a stripe not to split $1- P_{split}$ till a time t, as a function of its time or age, showing a non-exponential tail
(b) Stripes at $Re = 710$, a little above $Re_{c}$, exhibiting cycles of growth and splitting. The diagonal line indicates the growth of a stripe in absence of splitting, while the color gradient gives the probability for it to split, which increases with the length of the stripe. The gray curves give the evolution of the length of stripes with time, with one particular stripe shown in black The horizontal and vertical arrows indicate splitting events which 'reset' a stripes length and age. These events are numbered for the stripe shown in black - see text for details. After such splitting events, the stripe grows till a further splitting and so on, establishing a cycle of growth and splitting.  Note that here only splittings of stripe segments larger than $20h$ are shown for clarity. }
\label{fig:figure4}
\end{figure}

In contrast to the deterministic and constant growth at the stripe's downstream tip, the dynamics at the tail (i.e. upstream tip) is stochastic. Stripe segments of varying size split off and detach from the tail at irregular intervals. Some examples are marked by the red circles in Fig. \ref{fig:figure3a} and \ref{fig:figure3b}. The split pieces typically dissipate. However longer pieces occasionally survive and form a new stripe upstream and parallel to the parent stripe (see Supplementary Fig.3 for an example).
To quantify this, we look a the probability of splitting as a function of the stripe length. The constant growth rate of a stripe in between splittings can be used to convert the length into a time, which may be interpreted as the intrinsic time or age of a stripe, where time t = 0 corresponds to the extrapolated time where the stripe length is zero. The splitting statistics are shown in Fig.4a, where the probability of
a stripe to stay intact (i.e not split) is shown as a function of length / age of the stripe.
 Stripes shorter than 200h have a splitting probability close to zero. However as time proceeds and their lengths continue to increase, stripes appear to loose their structural stability and the probability to split increases, resulting in a non-exponential tail in Fig 4a. Thus, also stripe splitting is not memoryless, but an aging process.

The coupling between the growth of a stripe and its splitting probability leads to a cyclic process depicted in Fig. 4b. The length evolution is shown for several randomly chosen stripes (grey curves) at Re = 710. A representative stripe is highlighted in black. Starting from point 1 it initially grows at a constant rate (diagonal line in Fig. 4b), driven by the addition of streaks at the downstream tip. With the stripe's length also its probability to split increases (shown by the color gradient) and in this case the splitting eventually occurs at  point 2. Splitting resets this stripe to a shorter length (point 3) which can be interpreted as a younger stripe age and the growth resumes until the next splitting occurs (point 4). The stripe is ‘re-set’ to point 5, and grows until point 6 (at which point it exited the test section.
Corresponding cyclic dynamics (light grey curves) are exhibited by all stripes regardless of initial conditions.
The interplay between stripe elongation and splitting, causes the overall proliferation of turbulence to be two dimensional.

Aging dominates stripe dynamics below as well as above the critical point and both aspects are closely connected to the deterministic growth at the stripes' downstream tip, an aspect that is crucially missing in many previous studies of channel flow.

Memoryless stochastic dynamics is a defining characteristic of shear flow turbulence, It is a signature of the underlying chaotic saddle and in a spatio-temporal context this Markov property allows to link the problem to directed percolation. Many recent studies have confirmed this scenario suggesting its robustness. On the other hand examples from astrophysical and cardiovascular flows have shown  that subtle geometrical changes and imperfections can entirely alter the flows' stability and the nature of the transition\cite{Ji2006,Xu2020}. In channel flow it is not imperfections but on the contrary the approach towards the classic problem's (i.e. unconstrained) geometry which challenges the robustness of the transition model.

\section{Experiments}
\subsection{Experimental Setup}
 The experiments are carried out in a large aspect ratio channel which consists of two plates separated by a narrow gap of $2h = 1 \pm 0.03$ mm (Fig.S\ref{fig:ChannelSetup}). The bottom plate is a 10 mm thick, finely polished aluminum plate, while the top is 10 mm thick float glass. The gap is maintained by 1mm thick steel strips, which are clamped between the glass and steel plates and form the side-walls of the set up. The size of the channel is $(Lx; Ly; Lz) = (4000h; 2h; 490h)$, where $L_x$, $L_y$ and $L_z$ are lengths in the streamwise, wall-normal and spanwise directions respectively.  The working fluid is water, supplied from a continually overflowing reservoir, which is kept at the height of 21 m above the channel. The water level in the reservoir is maintained within $\pm 2$ cm, ensuring the pressure-head that drives the flow to be constant within $\pm 0.1\%$. The excess pressure drop due to an isolated stripe  is less than 0.2$\%$ of the total pressure head and hence the flow rate remains constant during stripe decays or splittings within a fraction of a percent.

\begin{figure}[h]
 \includegraphics[width=0.99\textwidth]{SI_figure1.eps}%
 \caption{The Channel set-up}
 \label{fig:ChannelSetup}
\end{figure}

Prior to entering the channel the water passes through a settling chamber including a porous barrier which breaks up large eddies. It then enters the channel proper via a convergence with an area ratio of 90:1. Unless perturbed, the flow in the channel remains laminar over the entire Reynolds number range investigated here.

The water temperature is measured just before the convergence by a precision resistance thermometer. The viscosity of the water $\nu$ is then determined using a fit to standard temperature-viscosity tables for water. The flow rate is measured by a magnetic flowmeter (ABB) installed in the main supply pipe leading to the channel, from which the bulk velocity $u_{bulk}$ is calculated.

To permit comparison with earlier studies, quantities are non-dimensionalised with the length scale h (half-gap) and 1.5 times the bulk velocity $u_{bulk}$  (which corresponds to the centre line velocity of the laminar flow at the same mass flow rate).

\subsection{Perturbation techniques}

To study the evolution of isolated turbulent stripes, an effective mechanism to generate them is needed. For $Re>900$, almost any standard perturbation mechanism (static obstacle, moving obstacle, jets injected through one or more holes in the channel wall) can efficiently trigger turbulence. Although for the lower Re range studied here ($600 < Re <900$), many of these do not work anymore, two methods were found to be efficient in generating isolated stripes. In one, perturbing the flow over a certain minimum span-wise ($\sim20h$) width successfully triggered stripe turbulence. In the present study we realized this by placing a ferromagnetic obstacle $100h$ downstream of the inlet which can be actuated by impulsively moving an externally placed magnet across the upper channel surface. This extended perturbation generates a seed which evolves into a stripe.

In the other method, a stripe is generated at a Re higher than the target Re,  where stripes are easily generated by a localized or extended perturbation (as described above). Next the Reynolds number was quenched to the target value. Thus was accomplished with a bypass loop in the supply pipe leading to the channel, which was normally closed by a solenoid valve. Prior to triggering the stripe, the flow-rate is adjusted to the target Re, with the bypass valve being closed.  The valve is then opened, reducing the resistance in the supply pipes, and hence increasing the flow-rate and hence the Re by around 100 - 150. After the stripe is triggered at this higher Re, it is allowed to develop for a time of around $300 h/u_{bulk}$. Subsequently the solenoid valve is closed, quenching the flow down to the target Re. Readings are taken after a further $250 h/u_{bulk}$ to allow the flow to adjust to the change in Re.
In comparison to the quench, the generation of stripes by an extended perturbation becomes increasingly inefficient for the lowest Re studied here ($Re <700$).

\subsection{Flow visualization}
The channel is illuminated with the help of LED lights installed parallel to the channel axis on both sides of the channel, approximately $1.4$ m  above it.
The flow was visualized using reflective flakes $10 - 40 \mu$ in size (Eckart SYMIC C001 reflective mica particles). These flakes tend to orient along the shear and hence allow to clearly distinguish turbulent from laminar regions.
The flow structures are monitored in the test section with the help of three 4 Mpixel 100Hz cameras placed at a distance around $1.5$ m above the channel. In order to allow time for the quench or the decay of initial transients, the test section is located some distance downstream of the channel entrance, and flows are recorded from $1500h$ downstream of the channel entrance till $500h$ before the channel exit. The combined view-field of the three cameras is $2000h \times 490h$.
To monitor the evolution of a stripe, images are captured by the cameras at a set frequency, with the three cameras being simultaneously triggered each time by a TTL signal.

\subsection{Image Processing}

\begin{figure}[h]
 \includegraphics[width=0.8\textwidth]{ImageProcessing.eps}%
 \caption{Image Processing (a) Raw image from the cameras stitched together showing an isolated stripe. (b) Background subtracted image. (c) Image after filtering and sharpening. (d) Binary image after thresholding}
 \label{fig:Images}
\end{figure}

The images of the three cameras are merged to capture the flow field in the entire test section (Fig.S\ref{fig:Images}a) and further processed by subtracting images of the laminar background flow (Fig.S\ref{fig:Images}b)
Images are next low-pass filtered using a Wiener filter with a 5 x 5 pixel mask. This is sharpened using unsharp masking in order to increase the contrast between the edges of the stripe and its surroundings (Fig.S\ref{fig:Images}c). This can be further used for extracting quantities of interest. For instance, in order to determine the area occupied by a stripe, a simple threshold is subsequently applied to create a binary image, shown in Fig.S\ref{fig:Images}d. The pixel count of the black area is then the area occupied by the turbulent stripe. For a reasonable range of parameters used in the image processing (e.g. the threshold for binarizing the image), there is only a minor variation in derived quantities such as the area of the turbulent stripe. More importantly such changes do not affect the value of the critical point or the lifetime distributions etc.

The growth rate of turbulence was estimated as follows. For each Re, around 1000 stripes are generated. Each stripe is considered at a definite time after the perturbation, its area determined and then an average found over all the stripes at that Re at the same time, thus yielding an ensemble average of the turbulent fraction at that time. Doing this for different times gives  the evolution of the turbulent fraction as a function of time, from which the growth rate is determined. In flows where the proliferation and decay of turbulence are memoryless, the typical time scales of these processes are well defined and the critical point is then determined as the Re where these timescales are in balance. By contrast, here, the critical point is determined as the Re where the growth rate changes from negative to positive.

\subsection{Stripe splitting}

 Small stripe segments are randomly shed from the stripe's upstream (trailing) edge. One such event each is shown by red circles in Fig 3(a) and (b) of the main paper. At low Re ($\lesssim 675$), such pieces typically decay. At larger Re however, some of these develop into an independent stripe parallel to, and upstream of the original parent stripe. An example of this splitting process is shown in Fig.S\ref{fig:Split}.

\begin{figure}[h]
 \includegraphics[width=0.99\textwidth]{Split.eps}%
 \caption{A sequence of three images at  $Re = 750$ showing stripe splitting. Consecutive images are 800 advective time units apart. A portion at the upstream tip (trailing edge) breaks off (b) and grows into an independent stripe upstream of the original one (c) }
 \label{fig:Split}
\end{figure}

\section{Numerical Simulations}
The numerical simulations in the present study are carried out using a modified version of the spectral code \textit{openpipeflow} \cite{openpipeflow}, which was adapted to simulate channel flow  \cite{xiao2020growth} in a rectangular box.

$X',Y',Z'$ are the streamwise, spanwise and wall-normal directions, with the associated unit vectors being $e_x',e_y',e_z'$ respectively. When scaled with half channel height $h$ and centre-line velocity $U_{cl}$ of the parabolic profile with the same mass flux, the laminar base flow can be expressed as  ${\bm U}=(1-y^2){\bm e}_{x'}$ , and  the Navier-Stokes equations for the fluctuating components  $(\bm{u},p)$  around base flow ${\bm U}$ take the form

\begin{align}
\frac{\partial \bm{u}}{\partial t} + (\bm{u} \cdot  \nabla) \bm{u} +(\bm{u} \cdot  \nabla) \bm{U} + (\bm{U} \cdot  \nabla) \bm{u} \nonumber \\
= -\nabla p + \frac{1}{Re} \nabla^2 \bm{U}  + \frac{1}{Re} \nabla^2 \bm{u} + {\bm f}(t)\\
\nabla \cdot \bm{u}= &0
\label{eq:NS_perturb}
\end{align}

 Here, the Reynolds number is defined as $Re=U_{cl}h/\nu$, where $\nu$ is the kinematic viscosity of the fluid.   $\bm{f}(t)=f(t)\bm{e}_{x'}$ is a time-dependent forcing term that represents an imposed streamwise pressure gradient. The mass flux in the streamwise direction $\bm{e}_{x'}$ is kept constant by changing the amplitude of the forcing term at every time-step.

No slip boundary conditions are imposed on the domain walls
\begin{align}
\bm u(x,\pm 1,z) &= 0 \nonumber \\
\label{eq:NoSlip_BC}
\end{align}

\subsection{Computational domain}
The computations are carried out in a rectangular box shaped domain which can be tilted with respect to the streamwise direction \cite{barkley2005computational,tuckerman2014turbulent} at an angle $\theta$, where $0^{\circ} \leq \theta \leq 90^{\circ}$.

The sides of the rectangular box  $x$ and $z$ have associated unit vectors ${\bm e}_{x}$ and ${\bm e}_{z}$ respectively. The usual streamwise, wall-normal and spanwise directions have following relation with the unit vectors associated with the domain directions.
\begin{align}
 \hat{{\bm e}}_{x'}&=  \cos{\theta}\hat{{\bm e}}_{x} + \sin{\theta}\hat{{\bm e}}_{z}  \\
 \hat{{\bm e}}_{z'}&= -\sin{\theta}\hat{{\bm e}}_{x} + \cos{\theta}\hat{{\bm e}}_{z} \\
 \hat{{\bm e}}_{y'}&= \hat{{\bm e}}_{y}
  \label{eq:tilt_coord}
\end{align}
 A constant mass flux is maintained in the streamwise direction and periodic boundary conditions are imposed on the faces normal to ${\bm e}_{x}$ and ${\bm e}_{z}$.

\begin{align}
\bm u(x,y,z)&= \bm u(x+L_{x},y,z) \nonumber \\
\bm u(x,y,z)&= \bm u(x,y,z+L_{z})
\label{eq:periodic_BC}
\end{align}
where, $L_x$ and $L_z$ are lengths of the domain in the $x$ and $z$ directions respectively.

The solver uses Fourier modes in the two periodic directions and  finite-difference in the wall-normal direction. The velocity field is decomposed as
\begin{align}
\bm u(x,y,z,t) = \sum_{k=-K}^{K} \sum_{m=-M}^{M} \hat{\bm u}_{k,m}(y,t)e^{i(\alpha k_{x}x + \beta m_{z} z)}
\label{eq:velDecompose1}
\end{align}
where $k$ and $m$ are Fourier modes in the $x$ and $z$ directions respectively, $\alpha= 2\pi/L_{x}$, $\beta= 2\pi/L_{z}$. $L_{x}$ ,$L_{z}$ are the lengths of the domain in $x$ and $z$ directions respectively. The time integration is performed using a second-order backward differentiation for linear terms and the Adam-Bashforth method for the nonlinear terms.

 The evolution of the stripes is monitored by defining a scalar observable - perturbation kinetic energy  $E(t)$ defined as

\begin{align}
E(t)=\frac{1}{L_{x}L_{y}L_{z}}\int_{-1}^{1}\left(\int_{0}^{L_{z}}\int_{0}^{L_{x}}|{\bm u}|^2dxdz\right)dy.
\end{align}

\subsection{Large domain with $\theta=0^{\circ}$}
A fully localized stripe can be simulated only in a large enough domain. We studied the lifetimes of the stripes at $Re=620$ in a square domain with no tilt i.e. $\theta=0^{\circ}$ with size $(L_{x},L_{y}, L_{z})=(400, 2, 400)$ and resolution $(N_{x},N_{y},N_{z})=(1536, 64, 1536)$.

\begin{figure}[h]
 \includegraphics[width=0.99\textwidth]{SI_figure2.eps}%
 \caption{Survival probability of fully localised stripes at $Re=620$. For different values of cutoffs used to determine lifetimes, the qualitative behaviour does not change.}
 \label{fig:Cutoffs}
\end{figure}

The resolution used is sufficient for the simulations at $Re=620$. The initial conditions for the lifetime study are uncorrelated snapshots from simulations at Re = 700 in the same domain. The lifetimes are determined by setting an appropriate cut-off on the perturbation kinetic energy. Once the perturbation KE falls below this threshold, the stripe is considered to have decayed. However, unlike the case of the partially localized stripes, there is no plateau in the energy before the decay. Rather the KE decays gradually (Fig.2b in the main paper). Hence, the lifetimes depend significantly on the threshold. The lifetimes for different threshold values are shown in Fig.S\ref{fig:Cutoffs}. Though the lifetimes do depend on the threshold, the qualitative behavior remains the same, exhibiting a non-exponential decay of survival probabilities.

\subsection{Tilted domain}
A partially localized stripe i.e. a stripe localized only in one direction and periodic in other direction is simulated in a rectangular domain tilted with respect to the streamwise direction at an angle of $\theta=45^{\circ}$. This tilt angle also corresponds to the stripe angle with respect to the streamwise direction \cite{paranjape2020oblique} and the angle of tilt chosen here agrees with the tilt angle of the stripes observed in the experiments and reported in \cite{tao2018extended}.

The domain size is $(L_{x},L_{y}, L_{z})=(30, 2, 100)$ with resolution $(N_{x},N_{y},N_{z})=(256, 49, 768)$.  The initial conditions for the lifetime study are taken from uncorrelated snapshots from simulations at $Re=600$. The stripe is considered to have decayed when the energy of fluctuations drops below a selected cutoff value e.g. $E=0.01$. Due to the rather clear plateau in energy before decay, the results do not strongly depend on the cut off chosen (see Fig.2 in the main paper).