# Binary patterns in the Prouhet-Thue-Morse sequence

J. Almeida  CMUP, Dep. Matemática, Faculdade de Ciências, Universidade do
Porto, Rua do Campo Alegre 687, 4169-007 Porto, Portugal<EMAIL_ADDRESS>and  O. Klíma  Dept. of Mathematics and Statistics, Masaryk University,
Kotlářská 2, 611 37 Brno, Czech Republic<EMAIL_ADDRESS>
###### Abstract.

We show that, with the exception of the words $a^{2}ba^{2}$ and $b^{2}ab^{2}$,
all (finite or infinite) binary patterns in the Prouhet-Thue-Morse sequence
can actually be found in that sequence as segments (up to exchange of letters
in the infinite case). We also identify the (finitely many) finite binary
patterns that appear non trivially, in the sense that they are obtained by
applying an endomorphism that does not map the set of all segments of the
sequence into itself.

###### Key words and phrases:

Prouhet-Thue-Morse sequence, pattern, infinite word, special word

###### 2010 Mathematics Subject Classification:

Primary: 68R15; Secondary: 11B85, 20M05, 20M35, 37B10

##  1\. Introduction

Let $\mu$ be the endomorphism of the free semigroup $\\{a,b\\}^{+}$ defined by
$\mu(a)=ab$ and $\mu(b)=ba$. Since $a$ is a prefix of $\mu(a)$, $\mu^{n}(a)$
is also a prefix of $\mu^{n+1}(a)$. Hence, the sequence $(\mu^{n}(a))_{n}$
determines a sequence of letters, or infinite word, whose prefix of length
$2^{n}$ is $\mu^{n}(a)$; we denote it $\mathbf{t}$. It is called the _Prouhet-
Thue-Morse sequence_ and it has been the object of extensive studies and
applications. It was first considered by Prouhet [12] in connection with a
problem in number theory, five decades later by Thue [20, 21] to exhibit
infinite words avoiding cubes and squares, and another two decades later by
Morse [11] as a discretized description of non-periodic recurrent geodesics in
surfaces of negative curvature. See [1] for a survey on this topic, including
several further connections with other branches of Mathematics. The first
author and other collaborators have previously studied the sequence
$\mathbf{t}$ in the framework of symbolic dynamics and its connections with
free profinite semigroups (see [3] and [4]). It was in fact an attempt to
construct a profinite semigroup with certain properties that prompted this
work, although no further references to profinite semigroups will be made in
this paper.

This paper concerns the study of binary patterns of $\mathbf{t}$, that is,
finite or infinite words $w$ over the alphabet $\\{a,b\\}$ for which there
exists an endomorphism $\varphi$ of the semigroup $\\{a,b\\}^{+}$ (naturally
extended to infinite words) such that the word $\varphi(w)$ can be found as a
block of consecutive letters of $\mathbf{t}$ (which we call a _segment_ of
$\mathbf{t}$). Since we need to identify concrete finite segments of
$\mathbf{t}$, a simple and efficient algorithm on how to compute them is
presented in Section 2.

Characterizations of binary patterns of $\mathbf{t}$ are due to Shur [16] and
D. Guaiana (unpublished work announced in [14, 15]). Our first main
contribution is a proof of the characterization attributed to D. Guaiana (but
also, independently, obtained by Shur in his thesis [18]) using Shur’s results
from [16]: with the exception of $a^{2}ba^{2}$ and $b^{2}ab^{2}$, the binary
patterns of $\mathbf{t}$ are its finite segments. Section 3 presents our proof
of this result.

The endomorphism $\mu$ and the endomorphism $\xi$ exchanging the letters $a$
and $b$ are easily seen to transform finite segments of $\mathbf{t}$ into
other such segments. In Section 4, we consider the problem of determining
which finite segments may only be transformed into other segments by
endomorphisms of $\\{a,b\\}^{+}$ that may be obtained by composition of $\mu$
and $\xi$. Such words are said to be _typical_ since we show that all but
finitely many finite segments of $\mathbf{t}$ are typical. We further
determine all atypical words. As an application of our results, we also
determine all infinite binary patterns of $\mathbf{t}$.

We conclude the paper with Section 5, where we propose the investigation of
the properties we established for $\mathbf{t}$ for arbitrary infinite words.

##  2\. Segments of $\mathbf{t}$

Note that $\\{ab,ba\\}$ is a code, in the sense that it generates a free
subsemigroup of $\\{a,b\\}^{+}$ and, therefore, $\mu$ is injective.

For an infinite word $w=a_{1}a_{2}\cdots$, by the _segments_ of $w$ we mean
the words of the form $a_{m}a_{m+1}\cdots a_{n}$ with $m\leqslant n$ or the
infinite words of the form $a_{m}a_{m+1}\cdots$. Note that, since
$\mu^{n+1}(a)=\mu^{n}(a)\mu^{n}(b)$, all factors of the words $\mu^{n}(b)$ are
segments of $\mathbf{t}$. It follows that a word $u\in\\{a,b\\}^{+}$ is a
segment of $\mathbf{t}$ if and only if so is the word that is obtained from
$u$ by interchanging the letters $a$ and $b$.

A word $w\in A^{+}$ is said to be _avoided_ by $\mathbf{t}$ if there is no
homomorphism $\varphi:A^{+}\to\\{a,b\\}^{+}$ such that $\varphi(w)$ is a
segment of $\mathbf{t}$. We also say that $w\in A^{+}$ is _unavoidable_ in
$\mathbf{t}$ if it is not avoided by $\mathbf{t}$; we then also say that $w$
is a _pattern_ of $\mathbf{t}$. For instance, it is well known that $a^{3}$
and $ababa$ are avoided by $\mathbf{t}$, which is also expressed by saying
that $\mathbf{t}$ is, respectively, cube-free and overlap-free [10].

The preceding notions are extended to infinite words by saying how
endomorphisms of $\\{a,b\\}^{+}$ are applied to infinite words. Given an
infinite word $w=a_{1}a_{2}\cdots$ over the alphabet $\\{a,b\\}$ and an
endomorphism $\varphi$ of $\\{a,b\\}^{+}$, we let $\varphi(w)$ be infinite
word obtained by concatenating the $\varphi(a_{i})$:
$\varphi(w)=\varphi(a_{1})\varphi(a_{2})\cdots$.

For a nonempty word $u$, let $t_{1}(u)$ denote its last letter.

The computation of the segments of $\mathbf{t}$ may be carried out easily in
view of the following proposition. The first part is an improved version of
[17, Corollary 1], although the same conclusion is in fact already established
in the proof of the cited statement. We present a proof for the sake of
completeness.

######  Proposition 2.1.

A word $w$ is a segment of $\mathbf{t}$ if and only if it is a factor of
$\mu^{n}(a)$, where $n=1$ if $|w|=1$, $n=3$ if $|w|=2$, and
$n=2+\lceil\log_{2}(|w|-1)\rceil$ otherwise. Moreover, for every integer
$k\geqslant 3$, the value $n=2+\lceil\log_{2}(k-1)\rceil$ is minimum for
$\mu^{n}(a)$ to admit as factors all segments of $\mathbf{t}$ of length $k$.

###### Proof.

Since $\mathbf{t}$ is cube-free, the cases where $|w|\leqslant 3$ are easily
verified by inspection. Suppose that $w$ is a segment of $\mathbf{t}$ which we
may assume to be of length at least 4. Then, $w$ is a factor of $\mu^{n}(a)$
for some positive integer $n$. Take $n$ to be minimum with that property. If
$m$ is the minimum positive integer such that $w$ is a factor of $\mu^{m}(x)$
for some letter $x$, then either only $x=b$ can play that role and $n=m+1$, or
$x=a$ may play it and $n=m$. We need to show, respectively, that $m\leqslant
1+\lceil\log_{2}(|w|-1)\rceil$ or $m\leqslant 2+\lceil\log_{2}(|w|-1)\rceil$.
Since $|w|\geqslant 4$, we may assume that $m\geqslant 4$ for, otherwise, the
inequality $m\leqslant 1+\lceil\log_{2}(|w|-1)\rceil$ holds trivially.

Let $\mu(x)=xy$. As $\mu^{m}(x)=\mu^{m-1}(xy)$ and $m$ is minimum, there must
be a nontrivial factorization $w=w_{1}w_{2}$ with $w_{1}$ a suffix of
$\mu^{m-1}(x)$ and $w_{2}$ a prefix of $\mu^{m-1}(y)$.

If one of the factors $w_{1}$ or $w_{2}$ has length greater than $2^{m-2}$,
then we must have $|w|>2^{m-2}+1$, which implies that $m\leqslant
1+\lceil\log_{2}(|w|-1)\rceil$ and fulfills our aim. Thus, we may assume that
both $w_{1}$ and $w_{2}$ have length at most $2^{m-2}$. Now, we have
$\mu^{m}(x)=\mu^{m-2}(xyyx)$ and $w$ is a factor of
$\mu^{m-2}(yy)=\mu^{m-3}(yxyx)$. If $w$ is a factor of either $\mu^{m-3}(xyx)$
or $\mu^{m-3}(yxy)$, then it is also a factor of
$\mu^{m-2}(xx)=\mu^{m-3}(xyxy)$ and, therefore, also of $\mu^{m}(y)$, so that
we are in the case $n=m$. On the other hand, by the minimality of $m$ the word
$w$ cannot be a factor of $\mu^{m-3}(xy)=\mu^{m-2}(x)$ or
$\mu^{m-3}(yx)=\mu^{m-2}(y)$, and so we have $|w|\geqslant 2^{m-3}+2$, which
implies that $n=m\leqslant 2+\lceil\log_{2}(|w|-1)\rceil$, as claimed. It
remains to consider the case where $w$ is a factor of neither $\mu^{m-3}(xyx)$
nor $\mu^{m-3}(yxy)$. Then there must be a factorization $w=s\mu^{m-3}(xy)z$
with $s$ and $z$ nontrivial words, so that $|w|\geqslant 2^{m-2}+2$, which
yields $m\leqslant 1+\lceil\log_{2}(|w|-1)\rceil$. This completes the proof of
the first part of the proposition.

To prove the last part of the proposition, first note that, for $k\geqslant
3$, the value of $f(k)=2+\lceil\log_{2}(k-1)\rceil$ is at least 3. We claim
that, for $n\geqslant 3$, there is a word of length $2^{n-3}+2$ that is a
segment of $\mathbf{t}$ but not a factor of $\mu^{n-1}(a)$. Noting that
$f(2^{n-3}+2)=n$, the result follows.

To establish the claim, consider the word
$w=t_{1}(\mu^{n-3}(b))\mu^{n-3}(a)b$. It is a segment of $\mathbf{t}$, in fact
a factor of
$\mu^{n}(a)=\mu^{n-2}(a)\mu^{n-3}(b)\cdot\mu^{n-3}(a)\cdot\mu^{n-1}(b)$ since
$b$ is the first letter of $\mu^{n-1}(b)$. It remains to show that $w$ is not
a factor of $\mu^{n-1}(a)=\mu^{n-3}(a)\mu^{n-3}(b)\mu^{n-3}(b)\mu^{n-3}(a)$.
Otherwise, since there are no overlaps in $\mathbf{t}$, $w$ must be a factor
of $\mu^{n-3}(bb)$. For $n=3$, this is clearly impossible since not even
$\mu^{n-3}(a)=a$ is a factor of $\mu^{n-3}(bb)=bb$. For $n>3$, we have

(1) $\displaystyle\mu^{n-3}(a)$ $\displaystyle=\mu^{n-4}(a)\mu^{n-4}(b)$ (2)
$\displaystyle\mu^{n-3}(bb)$
$\displaystyle=\mu^{n-4}(b)\mu^{n-4}(a)\mu^{n-4}(b)\mu^{n-4}(a).$

Since there are no overlaps in $\mu^{n-3}(bb)$, the only place where
$\mu^{n-3}(a)$ is found as a factor of $\mu^{n-3}(bb)$ is as the product of
the two middle factors in the factorization given in (2). Hence,
$w=t_{1}(\mu^{n-3}(b))\mu^{n-3}(a)b$ is not a factor of $\mu^{n-3}(bb)$ since,
for instance, $b$ is not the first letter of $\mu^{n-4}(a)$. ∎

For example, the segments of lengths 4 and 5 of $\mathbf{t}$ are the factors
of those lengths of $\mu^{4}(a)=abbabaabbaababba$. But, for instance, $aabb$
is a segment of $\mathbf{t}$ but not a factor of $\mu^{3}(a)=abbabaab$; it is
precisely the segment considered in the last part of the proof of Proposition
2.1.

Throughout the remainder of the paper, when we need to check whether a
concrete finite word is a segment of $\mathbf{t}$, without any further
reference we simply apply the algorithm given by Proposition 2.1, which is
linear in the length of the given word. We proceed similarly when we need to
compute all the segments of $\mathbf{t}$ of a given length.

Since every word of length $2^{n+2}-1$ must contain the factors $\mu^{n}(a)$
and $\mu^{n}(b)$, it follows from (the proof of) Proposition 2.1 that, for
$k\geqslant 3$, every segment of length $k$ of $\mathbf{t}$ is a factor of
every segment of $\mathbf{t}$ of length $\ell\leqslant
2^{2+1+\lceil\log_{2}(k-1)\rceil}-1\leqslant 16k-17$. The existence of such an
$\ell$ is the property known as _uniform recurrence_ of $\mathbf{t}$ and holds
for every sequence generated by iterating a primitive endomorphism of a free
semigroup [13, Proposition 5.2]. In the case of the Prouhet-Thue-Morse
sequence, the optimum value of $\ell$ is presented in [2, Example 10.9.3]: for
$k\geqslant 3$, we have $\ell=9\cdot 2^{r}+k-1$, where $r$ is the integer
determined by the inequalities $2^{r}+2\leqslant k\leqslant 2^{r+1}+1$. Note
that using the first inequality determining $r$, one gets the upper bound
$\ell\leqslant 10k-19$, which is better than our rough upper bound
$\ell\leqslant 16k-17$.

##  3\. Finite binary patterns

The following result plays a key role below.

######  Theorem 3.1 (Shur [16]).

The set of words of $\\{a,b\\}^{+}$ that are avoided by $\mathbf{t}$ is the
fully invariant ideal generated by the set

$\displaystyle\\{a^{3},ababa,a^{2}ba^{2}b,ab^{2}ab^{2},t_{1}(\mu^{k}(a))\mu^{k}(aba)a\
(k\geqslant 1),$ $\displaystyle\quad t_{1}(\mu^{m}(a))\mu^{m}(bab)a\
(m\geqslant 2)\\}.$

Moreover, the above is a minimal generating set for the fully invariant ideal
of the words avoided by $\mathbf{t}$.

The generators corresponding to $k=1$ and $k=2$ are, respectively,
$bab^{2}a^{2}ba$ and $a\,ab^{2}a\,ba^{2}b\,ab^{2}a\,a$; the generator
corresponding to $m=2$ is $a\,ba^{2}b\,ab^{2}a\,ba^{2}b\,a$ while, for $m=1$,
the word $t_{1}(\mu^{m}(a))\mu^{m}(bab)a=b^{2}a^{2}b^{2}a^{2}$ is avoided by
$\mathbf{t}$ but may be obtained for instance from the generator
$a^{2}ba^{2}b$ by mapping $a$ to $b$ and $b$ to $a^{2}$.

Another useful ingredient in our arguments is the following “synchronization”
result.

######  Lemma 3.2 ([8, Lemma 3.9]).

Let $X=\\{ab,ba\\}$ and consider $s\in X^{+}$ with $|s|\geqslant 4$. If $u$
and $v$ are words such that $usv\in X^{+}$ and $|u|$ is odd then $usv$ has an
overlap.

Since $\mathbf{t}$ has no overlaps, we conclude that $\mathbf{t}=usv$, with
$s$ as in Lemma 3.2 and $u$ a finite word, then $u$ has even length.

######  Corollary 3.3.

If there is a factorization $\mathbf{t}=u\mu^{n+1}(x)v$ where
$u\in\\{a,b\\}^{*}$, $x\in\\{a,b\\}$, and $n\geqslant 0$, then $u=\mu^{n}(w)$
and $v=\mu^{n}(z)$ for some word $w$ and infinite word $z$.

###### Proof.

We proceed by induction on $n$, the case $n=0$ being trivial as $\mu^{0}$ is
interpreted to be the identity function. Suppose that $n\geqslant 1$ and, by
symmetry, assume that $x=a$. Since
$\mu^{n+1}(a)=\mu^{n}(a)\mu^{n}(b)=\mu^{n-1}(abba)$, by the induction
hypothesis we know that $u=\mu^{n-1}(w)$ and $v=\mu^{n-1}(z)$, for some finite
word $w$ and infinite word $z$, and so $\mathbf{t}=wabbaz$. Lemma 3.2 then
implies that $w$ and $z$ belong to the image of $\mu$. Hence, $u$ and $v$
belong to the image of $\mu^{n}$. ∎

We say that a segment $u$ of $\mathbf{t}$ is _special_ if both $ua$ and $ub$
are segments of $\mathbf{t}$. The special segments of $\mathbf{t}$ have been
investigated by de Luca and Varricchio [8] with the purpose of counting the
number of segments of each given length. For our purposes, it suffices to
observe the following much simpler result.

######  Lemma 3.4 ([8, Lemma 3.6]).

If the word $w$ is special, then so is $\mu(w)$.

We say that two words are _suffix comparable_ if at least one of them is a
suffix of the other. The following lemma is the core key of our arguments.

######  Lemma 3.5.

Suppose that $u$ is a finite word such that $ua$ is unavoidable in
$\mathbf{t}$ and $ub$ is a segment of $\mathbf{t}$ but $ua$ is not. If
$n\geqslant 2$ and $u$ is suffix comparable with $\mu^{n}(a)$, then $u$ is
also suffix comparable with $\mu^{n+1}(b)$.

###### Proof.

Since $\mu^{n}(a)$ is a suffix of $\mu^{n+1}(b)=\mu^{n}(b)\mu^{n}(a)$, we may
assume that $\mu^{n}(a)$ is a suffix of $u$, say $u=u_{1}\mu^{n}(a)$. Consider
a concrete occurrence of $ub$ in $\mathbf{t}$: $\mathbf{t}=u_{0}ubv$, where
$u_{0}\in\\{a,b\\}^{*}$. Since $\mathbf{t}$ is recurrent, we may assume that
$|u_{0}|\geqslant 2^{n}$. By Corollary 3.3, we know that
$u_{0}u_{1}=\mu^{n-1}(u^{\prime})$ and $bv=\mu^{n-1}(v^{\prime})$ for some
word $u^{\prime}$ and infinite word $v^{\prime}$ (cf. Figure 1). Since
$\mu^{n-1}(v^{\prime})$ starts with $b$, so does $v^{\prime}$.

$u_{0}$$u$$bv$$u_{1}$$\mu^{n}(a)$$\mu^{n-1}(v^{\prime})$$\mu^{n-1}(u^{\prime})$
Figure 1. Some sgements of $\mathbf{t}$

Suppose first that $u^{\prime}$ ends with the letter $b$ and $|u|>3\cdot
2^{n-1}$. If $u^{\prime}$ ends in $ab$ then the word
$t_{1}(\mu^{n-1}(a))\mu^{n-1}(bab)a$ is a suffix of $ua$ which, in view of
Theorem 3.1, contradicts the assumption that $ua$ is unavoidable in
$\mathbf{t}$. On the other hand, if $u^{\prime}$ ends with $b^{2}$ then,
taking into account that $bv$ starts with $\mu^{n-1}(b)$, we conclude that
$\mathbf{t}=\mu^{n-1}(u^{\prime\prime}b^{2}ab^{2}v^{\prime\prime})$ for some
finite word $u^{\prime\prime}$ and infinite word $v^{\prime\prime}$. Since
$\mathbf{t}$ is a fixed point of the injective endomorphism $\mu$, it follows
that $b^{2}ab^{2}$ is a segment of $\mathbf{t}$, which we know is not the
case.

If $|u|\leqslant 3\cdot 2^{n-1}$, then $u$ is suffix of $\mu^{n-1}(xab)$ for a
letter $x$. By Lemma 3.4, as it is easy to check that $xab$ is special, so is
$\mu^{n-1}(xab)$. Hence, $u$ is special, contradicting the assumption that
$ua$ is not a segment of $\mathbf{t}$.

Thus, $u^{\prime}$ must end with $ba$, so that $u_{0}u$ ends with
$\mu^{n+1}(b)$. Since both $u$ and $\mu^{n+1}(b)$ are suffixes of $u_{0}u$,
they must be suffix comparable, thereby concluding the proof of the lemma. ∎

Similarly, one can prove the following lemma.

######  Lemma 3.6.

Suppose that $u$ is a finite word such that $ua$ is unavoidable in
$\mathbf{t}$ and $ub$ is a segment of $\mathbf{t}$ but $ua$ is not. If
$n\geqslant 2$ and $u$ is suffix comparable with $\mu^{n}(b)$ then $u$ is also
suffix comparable with $\mu^{n+1}(a)$.

Shur also observed in [16] that the word $a^{2}ba^{2}$ (and, therefore, also
$b^{2}ab^{2}$) is unavoidable in $\mathbf{t}$ but it is not a segment of
$\mathbf{t}$. Our first main result is that there are no other such examples,
thus providing an alternative characterization of two-letter words unavoidable
in t.

According to [15, Theorem 2] and [14, Theorem 3], the following theorem, which
is considered to be surprising, was first proved by D. Guaiana in 1996 but,
through private communication with A. Restivo, we learned that the proof was
never published and the manuscript appears to be lost. On the other hand, we
later learned from A. M. Shur that the next theorem also appears in his Ph.D.
thesis [18], dated 1997, which has never been published other than as a
document in the Russian State Library. Moreover, Shur observed that the result
can also easily be drawn from Theorem 3.1 using a characterization of the
finite words on the alphabet $\\{a,b\\}$ that are not segments of
$\mathbf{t}$, which is given in [19, Statement 1], a paper also in Russian.
Since all the proofs seem to be either lost or somewhat inaccessible in the
Russian literature, again we present our own proof for the sake of
completeness.

######  Theorem 3.7.

A word $w\in\\{a,b\\}^{+}$ is unavoidable in $\mathbf{t}$ if and only if it is
one of the words $a^{2}ba^{2}$ and $b^{2}ab^{2}$, or it is a segment of
$\mathbf{t}$.

###### Proof.

We proceed by induction on the length of $w$. In view of Theorem 3.1, it is
easy to check that the theorem holds for words of length at most 5. Assuming
inductively that the result holds for words of length $n$, let $w$ be a word
of length $n+1\geqslant 6$ that is unavoidable in $\mathbf{t}$. Since
interchanging the letters $a$ and $b$ does not affect either of the properties
of being unavoidable in $\mathbf{t}$ and being a segment of $\mathbf{t}$, we
may as well assume that $a$ is the last letter of $w$.

Let $w=ua$. Since $w$ is unavoidable in $\mathbf{t}$, so is $u$. Hence, by
induction hypothesis, $u$ may be found somewhere as a segment of $\mathbf{t}$.
Take such an occurrence of $u$ in $\mathbf{t}$ and let $x$ be the letter
immediately after it. We wish to show that there is such an occurrence of $u$
in $\mathbf{t}$ with $x=a$. Aiming at a contradiction, we may assume that
there is no such occurrence, that is, we always have $x=b$. Since $\mathbf{t}$
is recurrent, the segment $ux$ may be found in $\mathbf{t}$ as far as we wish,
so that we may continue prolonging it on the left as much as may be
convenient. Thus, we are assuming that $ua$ is unavoidable in $\mathbf{t}$ and
that $ub$ is a segment of $\mathbf{t}$ as long as desired but $ua$ is not a
segment of $\mathbf{t}$. However, we have to be careful because there is in
principle no assurance that such an extension of $u$ to the left retains the
property that $ua$ is unavoidable in $\mathbf{t}$.

Since $a^{3}$ is avoided by $\mathbf{t}$ and we are assuming that $x=b$, $u$
cannot end with $b^{2}$. We distinguish several cases according to the
termination of the word $u$.

If $u$ ends with $b$ then, by the above, it ends with $ab$. Suppose, more
precisely, that $u$ ends with $bab$. Since $w$ ends with $baba$ and $ababa$ is
avoided by $\mathbf{t}$, in fact $u$ must end with $b^{2}ab$. This situation
is impossible since we know that the suffix $b^{2}ab^{2}$ of $ub$ is not a
segment of $\mathbf{t}$.

Alternatively, assuming that $u$ ends with $b$, it must end with
$ba^{2}b=\mu^{2}(b)$. We may then apply successively Lemmas 3.5 and 3.6 to
deduce that there is some $n\geqslant 2$ such that $u$ is a suffix of
$\mu^{n}(a)$. By Lemma 3.4, it follows that $u$ is special, which contradicts
the assumption that $ua$ is not a segment of $\mathbf{t}$.

The next case we consider is that where $u$ ends with $aba$. Note that $u$
cannot end with $a^{2}ba$ for, otherwise, $w=ua$ ends with $ba^{2}ba^{2}$ and,
therefore, it cannot be unavoidable in $\mathbf{t}$. Also, $u$ cannot end with
$baba$ since $babab$ is not a segment of $\mathbf{t}$. Hence, $aba$ is not a
suffix of $u$.

Thus, assuming that $u$ ends with $ba$, it must end with $ab^{2}a=\mu^{2}(a)$.
We are then again led to a contradiction as above using Lemmas 3.5, 3.6 and
3.4. ∎

##  4\. Typical finite binary patterns

Recall that the endomorphism of $\\{a,b\\}^{+}$ switching the letters $a$ and
$b$ is denoted $\xi$. Since the finite segments of $\mathbf{t}$ are the finite
words that are factors of $\mu^{n}(a)$ for all sufficiently large $n$ and
$\mu^{n}(a)$ is a factor of $\mu^{n+1}(b)$, we may replace $a$ by $b$ in that
characterization of the segments of $\mathbf{t}$. Moreover, as $\xi$ commutes
with $\mu$, we conclude that the set of finite segments of $\mathbf{t}$ is
closed under applying the substitutions $\xi$ and $\mu$. By induction on $n$,
the words $\mu^{2n}(a)$ are palindromic, in the sense that they coincide with
the words read in the reverse order; this entails the well known fact that the
set of segments of $\mathbf{t}$ is closed under reversal.

We say that a word $w\in\\{a,b\\}^{+}$ is _atypical_ if it is a segment of
$\mathbf{t}$ and there is an endomorphism $\varphi$ of $\\{a,b\\}^{+}$ such
that $\varphi(w)$ is also a segment of $\mathbf{t}$ and $\varphi$ is not of
one of the forms $\mu^{n}$ or $\xi\circ\mu^{n}$ with $n\geqslant 0$. Segments
of $\mathbf{t}$ that are not atypical are said to be _typical_.

We say that a word is a _variant_ of another word $w$ if it may be obtained
from $w$ by applying reversal or $\xi$ or both. Note that the set of atypical
words is closed under taking factors and, by the above discussion, it is also
closed under taking variants.

The following result appears explicitly as [6, Proposition 3.3] but may
already be extracted from Thue’s work [21] (see [5, Chapter 3, Proposition
2.13]).

######  Proposition 4.1.

If $u^{2}$ is a segment of $\mathbf{t}$ then $u$ is one of the words
$\mu^{n}(a)$, $\mu^{n}(b)$, $\mu^{n}(aba)$ or $\mu^{n}(bab)$ for some
$n\geqslant 0$.

Yet another property of the Prouhet-Thue-Morse infinite word is the following
result which explains the above terminology.

######  Theorem 4.2.

Let $w\in\\{a,b\\}^{+}$ be a segment of $\mathbf{t}$ containing at least one
of the segments $aba$ and $bab$ along with all other segments of $\mathbf{t}$
of length 3. Then $w$ is typical.

###### Proof.

Suppose $\varphi$ is an endomorphism of $\\{a,b\\}^{+}$ such that $\varphi(w)$
is a segment of $\mathbf{t}$. By Proposition 4.1, since $\varphi(a^{2})$ and
$\varphi(b^{2})$ are square segments of $\mathbf{t}$, each of the words
$\varphi(a)$ and $\varphi(b)$ must be obtained by applying a power of $\mu$ to
one of the words $a,b,aba,bab$. Let then $\varphi(a)=\mu^{k}(u)$ and
$\varphi(b)=\mu^{\ell}(v)$, where $u,v\in\\{a,b,aba,bab\\}$. We may assume
that $k\leqslant\ell$ since, otherwise, we would consider the pair
$(\xi(w),\varphi\circ\xi)$ instead of $(w,\varphi)$. Then, we have the
factorization $\varphi=\mu^{k}\circ\psi$, where $\psi$ is endomorphism of
$\\{a,b\\}^{+}$ defined by $\psi(a)=u$ and $\psi(b)=\mu^{\ell-k}(v)$. Since
$\mu$ is injective and $\mathbf{t}$ is a fixed point of $\mu$, from the fact
that $\varphi(w)$ is a segment of $\mathbf{t}$, we conclude that so is
$\psi(w)$. One the other hand, since $\xi$ and $\mu$ commute, if $\psi$ is a
product of $\mu$ and $\xi$ then so is $\varphi$. Hence, we may assume that
$k=0$.

The mapping $\xi\circ\varphi$ has also the property that
$(\xi\circ\varphi)(w)$ is a segment of $\mathbf{t}$. Since
$(\xi\circ\varphi)(a)=\xi(\mu^{k}(u))=\mu^{k}(\xi(u))$, we may further assume
that $u$ is one of the words $a$ or $aba$. Since $aab$ is a factor of $w$ and
neither $a^{3}$ nor $a^{2}ba^{2}$ is a factor of $\mathbf{t}$, then $v$ cannot
start with $a$ and, therefore, it must be either $b$ or $bab$.

Consider first the case where $u=a$. Since $baa$ is a factor of $w$ but
$a^{3}$ is not a factor of $\mathbf{t}$, $\varphi(b)$ cannot end in $a$, which
implies that $\ell$ is even. If $\ell\geqslant 2$, then $\varphi(b)$ starts
with $\mu^{2}(b)=ba^{2}b$. But, since $a^{2}b$ is a factor of $w$, this
implies $a^{2}ba^{2}$ is a factor of $\varphi(w)$ and, therefore, a segment of
$\mathbf{t}$, which we know is not the case. Hence, we must have $\ell=0$. It
remains to rule out the case $\varphi(b)=bab$, which results from noting that
in that case, from the assumption that either $aba$ or $bab$ is a factor of
$w$, it follows that either $ababa$ or $bababab$ is a factor of $\varphi(w)$
while we know that $ababa$ is not a segment of $\mathbf{t}$.

Next, consider the case where $u=aba$. Since $a^{2}b$ is a factor of $w$ and
$ababa$ is not a segment of $\mathbf{t}$, $\varphi(b)$ cannot start with $ba$.
As $\mu^{\ell}(b)$ is a prefix of $\varphi(b)$, it follows that $\ell=0$ and
$\varphi(b)=b$. This leads to a similar situation as that considered at the
end of the preceding paragraph, with the letters $a$ and $b$ interchanged
which is, therefore, excluded. This completes the proof of the theorem. ∎

The assumption of Theorem 4.2 that a segment of $\mathbf{t}$ contains as
factors at least one of the words $aba$ and $bab$ along with all other
segments of length 3 of $\mathbf{t}$ holds for all segments of $\mathbf{t}$ of
length $10$, as may be easily checked by examining all segments of that
length. Hence, by Theorem 4.2 there are only finitely many atypical words.
Since there are 5 different segments of length 3 of $\mathbf{t}$ that are
supposed to appear in the word, no word with length shorter than 7 satisfies
the criterion of Theorem 4.2 and $ab^{2}a^{2}ba$ is a word of length 7 that
does satisfy it. On the other hand, the segment $a^{2}bab^{2}aba$, of length
9, fails to have the segment $ba^{2}$ as a factor.

The following result completes the above observations by giving the full
identification of atypical words.

######  Theorem 4.3.

Up to taking variants, the atypical words are the factors of the words
$aabab$, $abaaba$, and $aabbaab$.

###### Proof.

To check that all relevant words have been duly considered, the reader may
wish to refer to the diagram in Figure 2 later in the paper, where all
atypical words are represented.

The following is the complete list of segments of $\mathbf{t}$ of length 5:

$\displaystyle aabab,aabba,abaab,ababb,abbaa,abbab,$ $\displaystyle
baaba,baabb,babaa,babba,bbaab,bbaba.$

Note that all these words are variants of factors of at least one of the three
words in the statement of the theorem. Hence, by showing that those three
words are atypical, we obtain that so are all words of length up to 5.

We next indicate for each of the words in the statement of the theorem an
endomorphism $\varphi$ of $\\{a,b\\}^{+}$ not of the forms $\mu^{n}$ and
$\xi\circ\mu^{n}$ that maps it to a segment of $\mathbf{t}$:

  * •

$aabab$: $\varphi(a)=a$, $\varphi(b)=b^{2}aba^{2}b$;

  * •

$abaaba$: $\varphi(a)=a$, $\varphi(b)=b^{2}$;

  * •

$aabbaab$: $\varphi(a)=a$, $\varphi(b)=bab$.

The verification of all these statements amounts to routine calculations.

Less straightforward is the verification that there are no other atypical
words. Note that a word is typical if it has a typical factor. Hence, also
excluding variants and words that satisfy the criterion of Theorem 4.2, we
obtain the following reduced list of words remaining to be treated:

(3) $aababb,aabbab,abbaabba,ababba.$

We proceed to show that each word $w$ in the list (3) is typical. For that
purpose, assume that $\varphi$ is an endomorphism of $\\{a,b\\}^{+}$ such that
$\varphi(w)$ is a segment of $\mathbf{t}$.

In the first three cases, since $\varphi(a^{2})$ and $\varphi(b^{2})$ are
factors of $\varphi(w)$, we may start the argument using Proposition 4.1 as in
the proof of Theorem 4.2, assuming that $\varphi(a)$ is either $a$ or $aba$.

Consider first the case $\varphi(a)=a$. Since in the three cases, $aab$ is a
factor of $w$ but $\mathbf{t}$ is cube-free, $\varphi(b)$ must start with $b$.
Therefore, we may assume that $\varphi(b)=bv$ with $v\in\\{a,b\\}^{+}$. In all
the three first cases, since $abb$ is a factor of $w$, we get that $abvbv$ is
a factor of $\varphi(w)$ and this would provide an overlap in $\mathbf{t}$ if
$v$ ends with $a$. Hence $v$ ends with $b$. Since $\varphi(b)$ is of the form
$\mu^{n}(x)$ for some $x\in\\{a,aba,b,bab\\}$, we conclude that either
$\varphi(b)$ starts with $baab$ or it is $bab$. The first case is excluded
since $aabaa$ is then a factor of $\varphi(aab)$, whence also of $\varphi(w)$,
while it is not a segment of $\mathbf{t}$. The case $\varphi(b)=bab$ is also
excluded if $w$ is either $aababb$ or $aabbab$ since it leads to the overlap
$babab$ in the factor $\varphi(bab)$ of $\varphi(w)$. In case $w=abbaabba$,
one can simply check directly that $\varphi(w)$ is not a segment of
$\mathbf{t}$.

Still treating for the moment only the first three of the words in the list
(3), suppose next that $\varphi(a)=aba$. Again, as $aab$ is a factor of $w$
and $aabaa$ cannot be a factor of $\varphi(w)$, $\varphi(b)$ must start with
$b$. If it ends with $a$, then $a\varphi(bb)$ would be an overlap in
$\varphi(w)$ since $abb$ is a factor of $w$. Hence, $\varphi(b)$ starts and
ends with $b$. This is impossible in case $w$ has the factor $bab$ since it
would lead to the overlap $babab$ in $\varphi(w)$. This excludes the cases
where $w$ is the first or the second word in the list (3). So, we have
$w=abbaabba$. Then $\varphi(w)$ is a square segment of $\mathbf{t}$. By
Proposition 4.1, $\varphi(abba)$ is one of the words $\mu^{n}(x)$ with
$x\in\\{a,aba,b,bab\\}$. Since $n\leqslant 1$ gives a word that is too short
to be $\varphi(abba)$, we must have $n\geqslant 2$, in which case a simple
calculation shows that $\mu^{n}(x)$ cannot start with $aba$. This ends the
verification that the first three words in the list (3) are typical.

It remains to consider the word $w=ababba$. Here, we have two square factors
of $\varphi(w)$, namely the squares of $\varphi(b)$ and $\varphi(ab)$. By
Proposition 4.1 we know that there are words $x,y\in\\{a,aba,b,bab\\}$ and non
negative integers $m,n$ such that $\varphi(b)=\mu^{m}(x)$ and
$\varphi(ab)=\mu^{n}(y)$. In case $m>n$, comparing the lengths of the word
$\varphi(ab)$ and its factor $\varphi(b)$, we obtain the inequality
$2^{n}|y|>2^{m}|x|$, so that $3\geqslant|y|>2^{m-n}|x|\geqslant 2^{m-n}$. It
follows that $|x|=1$, $|y|=3$ and $m=n+1$. From the equalities
$\mu^{n}(y)=\varphi(ab)=\varphi(a)\mu^{n+1}(x)$ we then deduce that
$\varphi(a)=\mu^{n}(y_{1})$, where $y_{1}$ is the first letter of $y$, and
$\mu(x)=xy_{1}$. Since $abb$ is a factor of $w$,
$\varphi(abb)=\mu^{n}(y_{1}xy_{1}xy_{1})$ is a segment of $\mathbf{t}$, which
contradicts $\mathbf{t}$ being overlap free. Thus, we must have $n\geqslant
m$. Then $\varphi(a)$ must be of the form $\mu^{m}(z)$ where $z$ is a prefix
of $\mu^{n-m}(y)$. It follows that, as in the proof of Theorem 4.3, we may
then assume that $m=0$ and that $\varphi(b)$ is either $b$ or $bab$. Consider
first the case where $\varphi(b)=b$. Since $abba$ is a factor of $w$ but
$b^{3}$ is not, the word $\varphi(a)$ must start and end with the letter $a$
and we may assume that it is not reduced to $a$. Since $\varphi(a)=z$, we
conclude that $\varphi(a)$ must start with $abba$. Since $bba$ is a factor of
$w$, this yields the factor $bbabb$ of $\varphi(w)$, which is not possible
since $\varphi(w)$ is a segment of $\mathbf{t}$. Finally, the case
$\varphi(b)=bab$ is excluded since $\varphi(ab)=\mu^{n}(y)$ cannot end with
$bab$. This concludes the proof of the theorem. ∎

To facilitate the visualization of the set of atypical words, we give a
semigroup theoretical formulation. Although we do not go deep into it, the
reader unfamiliar with semigroup theory may prefer to skip these
considerations or refer to a standard textbook in the area such as [7, 9].

Let $S$ be the set of atypical words. We may define a multiplication on the
set $S^{0}=S\cup\\{0\\}$ as follows: for $u,v\in S$, $u\cdot v$ is $uv$ if
$uv$ is atypical and 0 otherwise; for all $s\in S^{0}$, $s\cdot 0=0\cdot s=0$.
Note that $S^{0}$ is the Rees quotient of $\\{a,b\\}^{+}$ by the ideal
consisting of the typical words together with the words that are not segments
of $\mathbf{t}$.

The diagram in Figure 2 represents $S^{0}$ as a partially ordered set for the
Green $\mathcal{J}$-order, in which an element $u$ lies above $v$ if and only
if $u$ is a factor of $v$. The words in bold are the lexicographic minima
among their variants; note that those that are atoms (which are underlined)
are precisely the words that were shown directly to be atypical in Theorem
4.3.

$\bm{a}$$\bm{aa}$$\bm{ab}$$b$$ba$$bb$$\bm{aab}$$baa$$\bm{aba}$$abb$$bab$$bba$$\bm{aaba}$$\bm{aabb}$$baab$$abaa$$\bm{abab}$$baba$$\bm{abba}$$babb$$bbaa$$bbab$$\underline{\bm{aabab}}$$baaba$$\bm{aabba}$$baabb$$\bm{abaab}$$babaa$$ababb$$abbaa$$abbab$$babba$$bbaab$$bbaba$$0$$\bm{aabbaa}$$baabba$$\underline{\bm{abaaba}}$$\bm{abbaab}$$babbab$$bbaabb$$\underline{\bm{aabbaab}}$$baabbaa$$abbaabb$$bbaabba$
Figure 2. The semigroup $S^{0}$

We conclude this section with another application of Theorem 4.2, this one
concerning infinite patterns of $\mathbf{t}$.

######  Corollary 4.4.

Let $w$ be an infinite word and suppose that there is an endomorphism
$\varphi$ of $\\{a,b\\}^{+}$ such that $\varphi(w)$ is a suffix of either
$\mathbf{t}$ or $\xi(\mathbf{t})$. Then $w$ is itself a suffix of either
$\mathbf{t}$ or $\xi(\mathbf{t})$.

###### Proof.

Since all segments of $w$ are unavoidable in $\mathbf{t}$ and they are all
extendable on the right, by Theorem 3.7 they are segments of $\mathbf{t}$.
Since the language of the segments of $\mathbf{t}$ defines a minimal subshift
[13, Proposition 5.2], it follows that $w$ and $\mathbf{t}$ have the same
segments. In particular, the word $a^{2}b^{2}a^{2}bab$ is a segment of $w$ and
it satisfies the assumption of Theorem 4.2. It follows that there is
$n\geqslant 0$ such that $\varphi=\mu^{n}$ or $\varphi=\xi\circ\mu^{n}$.
Again, since $\mu$ is injective and both $\mathbf{t}$ and $\xi(\mathbf{t})$
are fixed by $\mu$, the result follows. ∎

The somewhat different formulation for finite and infinite segments (compare
Theorem 3.7 with Corollary 4.4) is fully justified by the following result.

######  Proposition 4.5.

The infinite words $\mathbf{t}$ and $\xi(\mathbf{t})$ have no common suffix.

###### Proof.

Suppose on the contrary that the infinite word $w$ is a suffix of both
$\mathbf{t}$ and $\xi(\mathbf{t})$. Since every suffix of $w$ is again a
suffix of both $\mathbf{t}$ and $\xi(\mathbf{t})$, we may assume that
$\xi(\mathbf{t})=\mu^{n}(b)w$. Since $\xi(\mathbf{t})$ is a fixed point of
$\mu$ and $\mu$ is injective, we conclude that there is some infinite word $x$
such that $w=\mu^{n}(x)$ and $\xi(\mathbf{t})=bx$. Note that $x=a\mu(y)$ for
an infinite word $y$. Let $u$ be a finite word such that
$\mathbf{t}=uw=u\mu^{n}(a)\mu^{n+1}(y)$. By Corollary 3.3, we know that
$u\mu^{n}(a)=\mu^{n}(v)$ for some word $v$. Hence,
$\mathbf{t}=\mu^{n}(v\mu(y))$ and so, as above, we must have
$\mathbf{t}=v\mu(y)$. Since $x$ starts with $aababba$, $\mu(y)$ starts with
$abab$ and so $v$ must end with $ba$. This shows that
$\xi(\mathbf{t})=ba\mu(y)$ must be a suffix of $\mathbf{t}$. Interchanging the
roles of $\mathbf{t}$ and $\xi(\mathbf{t})$, we similarly conclude that
$\mathbf{t}$ is a suffix of $\xi(\mathbf{t})$. Hence, $\mathbf{t}$ is a proper
suffix of itself, which implies that it is periodic; this contradicts, for
instance, the fact that $\mathbf{t}$ is cube-free. ∎

##  5\. Final remarks and problems

For an infinite word $\mathbf{w}$ over a finite alphabet $A$, let
$L(\mathbf{w})$ be the language consisting of its finite segments. Note that
the automorphisms of the semigroup $A^{+}$ permute the letters of $A$; we call
them _letter exchanges_. The language obtained from $L(\mathbf{w})$ by
applying all possible letter exchanges is denoted $\bar{L}(\mathbf{w})$. We
denote by $\mathcal{E}(\mathbf{w})$ the set of all endomorphisms $\varphi$ of
$A^{+}$ such that $\varphi(L(\mathbf{w}))\subseteq L(\mathbf{w})$. The set
$\bar{\mathcal{E}}(\mathbf{w})$ is similarly defined using
$\bar{L}(\mathbf{w})$ instead of $L(\mathbf{w})$. Note that
$\mathcal{E}(\mathbf{w})$ is a submonoid of $\bar{\mathcal{E}}(\mathbf{w})$
which in turn is a submonoid of the monoid $\mathop{\mathrm{End}}(A^{+})$ of
all endomorphisms of the semigroup $A^{+}$.

The following is an immediate consequence of Theorem 4.2.

######  Corollary 5.1.

The monoid $\mathcal{E}(\mathbf{t})=\bar{\mathcal{E}}(\mathbf{t})$ is
generated by the set $\\{\xi,\mu\\}$. In particular, it is finitely
generated.∎

Corollary 5.1 is intimately related with a result of Thue (see [5, Chapter 3,
Theorem 2.16]) that characterizes the set of the so-called _overlap-free
morphisms_ , that is, endomorphisms of $\\{a,b\\}^{+}$ that map the set of all
overlap-free words into itself, namely as the monoid generated by
$\\{\xi,\mu\\}$. In fact, in view of another result of Thue (see [5, Chapter
3, Theorem 2.15]), all (overlap-free) words that can be arbitrarily prolonged
in both directions to overlap-free words are segments of $\mathbf{t}$. It
follows that overlap-free morphisms belong to $\mathcal{E}(\mathbf{t})$ and so
Corollary 5.1 immediately yields Thue’s necessary condition for overlap-free
morphisms. That the condition is also sufficient is given by another result of
Thue (see [5, Chapter 3, Lemma 2.2]). It does not appear to be immediately
obvious how to deduce Corollary 5.1 from Thue’s results.

Theorems 3.7 and 4.2, together with Corollary 5.1 may be regarded as three
finiteness properties of the Prouhet-Thue-Morse sequence. It is natural to ask
which infinite words possess such finiteness properties. More precisely, we
propose the following problems.

######  Problem 1.

Which infinite words $\mathbf{w}$ have the property that, up to finitely many
exceptions, the patterns of $\mathbf{w}$ on the same alphabet are obtained
from its segments up to an exchange of letters?

######  Problem 2.

For which infinite words $\mathbf{w}$ is the monoid $\mathcal{E}(\mathbf{w})$
finitely generated? Similar question for $\bar{\mathcal{E}}(\mathbf{w})$.

We say that a finite segment $u$ of $\mathbf{w}$ is _$\mathbf{w}$ -atypical_
if there is some endomorphism $\varphi\notin\bar{\mathcal{E}}(\mathbf{w})$ of
$A^{+}$ such that $\varphi(u)$ is also a segment of $\mathbf{w}$.

######  Problem 3.

Which infinite words $\mathbf{w}$ have only finitely many
$\mathbf{w}$-atypical segments?

A negative example for Problem 1 is provided by the _Fibonacci infinite word_
, which is the only fixed point $\mathbf{f}$ of the endomorphism _$\phi$_ of
$\\{a,b\\}^{+}$ defined by $\phi(a)=ab$ and $\phi(b)=a$. That there are
infinitely many finite binary patterns of $\mathbf{f}$ that are not segments
of $\mathbf{f}$ was proved in [15] (see also [14]), where it is also shown
that there are Sturmian infinite words that admit as patterns all segments of
all Sturmian infinite words. Recall that an infinite word is _Sturmian_ if it
has exactly $n+1$ segments of each length $n\geqslant 1$. We do not know
whether $\mathcal{E}(\mathbf{f})$ is generated by $\varphi$ and
$\bar{\mathcal{E}}(\mathbf{f})$ is generated by $\varphi$ and $\xi$. We also
do not know whether the set of $\mathbf{f}$-atypical words is finite.

Problem 1 was raised in [15] for binary infinite words that are either fixed
points of endomorphisms or of linear complexity. In the same paper, it is
observed that if $\mathbf{w}$ is an infinite word with all elements of $A^{+}$
as segments (which may be obtained for instance by concatenating all the words
in a sequence enumerating the elements of $A^{+}$), then obviously
$\mathbf{w}$ is a positive example for Problem 1. Note that
$\mathcal{E}(\mathbf{w})=\bar{\mathcal{E}}(\mathbf{w})=\mathop{\mathrm{End}}(A^{+})$
and it is easy to see that $\mathop{\mathrm{End}}(A^{+})$ is not finitely
generated: for the endomorphisms that maps each letter to itself, except for
one letter $a$ that is mapped to $a^{p}$, where $p$ is prime, the only
elements of $\mathop{\mathrm{End}}(\mathbf{w})$ that are factors of it are the
letter exchanges and the factors of which it is also a factor. From the
preceding observation it also follows that there are no $\mathbf{w}$-atypical
words. Thus, $\mathbf{w}$ is a negative example for Problem 2 and a positive
example for Problem 3.

## Acknowledgments

The first author acknowledges partial support by CMUP (UID/MAT/ 00144/2019)
which is funded by FCT (Portugal) with national (MATT’S) and European
structural funds (FEDER) under the partnership agreement PT2020. The work was
carried out at Masaryk University, whose hospitality is gratefully
acknowledged, with the support of the FCT sabbatical scholarship
SFRH/BSAB/142872/2018. The second author was supported by Grant 19-12790S of
the Grant Agency of the Czech Republic.

## References

  * [1] J.-P. Allouche and J. Shallit, _The ubiquitous Prouhet-Thue-Morse sequence_ , Sequences and their applications (Singapore, 1998), Springer Ser. Discrete Math. Theor. Comput. Sci., Springer, London, 1999, pp. 1–16.
  * [2] by same author, _Automatic sequences: Theory, applications, generalizations_ , Cambridge University Press, 2003.
  * [3] J. Almeida and A. Costa, _Presentations of Schützenberger groups of minimal subshifts_ , Israel J. Math. 196 (2013), 1–31.
  * [4] J. Almeida, A. Costa, R. Kyriakoglou, and D. Perrin, _Profinite semigroups and symbolic dynamics_ , 2019, A book in preparation.
  * [5] J. Berstel, _Axel Thue’s papers on repetitions in words: a translation_ , Tech. report, Université du Québec à Montréal, 1995, Publications du LaCIM 20, available at http://www-igm.univ-mlv.fr/~berstel/Articles/1994ThueTranslation.pdf.
  * [6] S. Brlek, _Enumeration of factors in the Thue-Morse word_ , Discrete Appl. Math. 24 (1989), 83–96.
  * [7] A. H. Clifford and G. B. Preston, _The algebraic theory of semigroups_ , vol. I, Amer. Math. Soc., Providence, R.I., 1961.
  * [8] A. de Luca and S. Varricchio, _Some combinatorial properties of the Thue-Morse sequence and a problem in semigroups_ , Theor. Comp. Sci. 63 (1989), no. 3, 333–348.
  * [9] J. M. Howie, _Fundamentals of semigroup theory_ , London Mathematical Society Monographs. New Series, vol. 12, The Clarendon Press, Oxford University Press, New York, 1995.
  * [10] M. Lothaire, _Combinatorics on words_ , Addison-Wesley, Reading, Mass., 1983\.
  * [11] H. M. Morse, _Recurrent geodesics on a surface of negative curvature_ , Trans. Amer. Math. Soc. 22 (1921), 84–100.
  * [12] E. Prouhet, _Mémoire sur quelques relations entre les puissances des nombres_ , C. R. Acad. Sci. Paris 33 (1851), 31.
  * [13] M. Queffélec, _Substitution dynamical systems—spectral analysis_ , second ed., Lect. Notes in Math., vol. 1294, Springer-Verlag, Berlin, 2010.
  * [14] A. Restivo and S. Salemi, _Binary patterns in infinite binary words_ , Formal and natural computing, Lect. Notes in Comput. Sci., vol. 2300, Springer, Berlin, 2002, pp. 107–116.
  * [15] by same author, _Words and patterns_ , Developments in language theory (Vienna, 2001), Lect. Notes in Comput. Sci., vol. 2295, Springer, Berlin, 2002, pp. 117–129.
  * [16] A. M. Shur, _Binary words avoided by the Thue-Morse sequence_ , Semigroup Forum 53 (1996), 212–219.
  * [17] by same author, _Overlap-free words and Thue-Morse sequences_ , Int. J. Algebra Comput. 6 (1996), 353–367.
  * [18] by same author, _Algebraic and combinatorial properties of rational languages_ , Ph.D. thesis, Ural State University, 1997, In Russian.
  * [19] by same author, _Combinatorial complexity of rational languages_ , Diskretn. Anal. Issled. Oper. 12 (2005), no. 2, 78–99, In Russian.
  * [20] A. Thue, _Über unendlichen zeichenreihen_ , Kra. Vidensk. Selsk. Skrifter, I. Mat. Nat. Kl. (1906), no. 7, 1–22.
  * [21] by same author, _Über die gegenseitige Lage gleicher Teile gewisser Zeichenreihen_ , Kra. Vidensk. Selsk. Skrifter, I. Mat. Nat. Kl. (1912), no. 1, 1–67.


# Parallel-in-Time Integration of the Landau-Lifshitz-Gilbert Equation with
the Parallel Full Approximation Scheme in Space and Time

Robert Kraft  Research Platform MMM Mathematics - Magnetism - Materials,
University of Vienna, Vienna, Austria  Vienna Doctoral School in Physics,
University of Vienna, Vienna, Austria  Sabri Koraltan  Physics of Functional
Materials, University of Vienna, Vienna, Austria  Vienna Doctoral School in
Physics, University of Vienna, Vienna, Austria  Markus Gattringer  Physics of
Functional Materials, University of Vienna, Vienna, Austria  Vienna Doctoral
School in Physics, University of Vienna, Vienna, Austria  Florian Bruckner
Research Platform MMM Mathematics - Magnetism - Materials, University of
Vienna, Vienna, Austria  Physics of Functional Materials, University of
Vienna, Vienna, Austria  Dieter Suess  Research Platform MMM Mathematics -
Magnetism - Materials, University of Vienna, Vienna, Austria  Physics of
Functional Materials, University of Vienna, Vienna, Austria  Claas Abert
Research Platform MMM Mathematics - Magnetism - Materials, University of
Vienna, Vienna, Austria  Physics of Functional Materials, University of
Vienna, Vienna, Austria

###### Abstract

Speeding up computationally expensive problems, such as numerical simulations
of large micromagnetic systems, requires efficient use of parallel computing
infrastructures. While parallelism across space is commonly exploited in
micromagnetics, this strategy performs poorly once a minimum number of degrees
of freedom per core is reached. We use magnum.pi, a finite-element
micromagnetic simulation software, to investigate the Parallel Full
Approximation Scheme in Space and Time (PFASST) as a space- and time-parallel
solver for the Landau-Lifshitz-Gilbert equation (LLG). Numerical experiments
show that PFASST enables efficient parallel-in-time integration of the LLG,
significantly improving the speedup gained from using a given number of cores
as well as allowing the code to scale beyond spatial limits.

††preprint: APS/123-QED

##  I Introduction

Computationally expensive simulations of large micromagnetic systems can be
significantly accelerated by carrying out computations in parallel.
Parallelism in space via distributed meshes or grids is commonly exploited in
micromagnetic simulation codes, for instance in CPU software such as OOMMF
[1], mumax [2] or GPU codes like magnum.np [3] or BORIS [4]. However, the
efficiency of this strategy drops sharply when the number of degrees of
freedom per core gets too small, as the communication overhead begins to
dominate the computational cost. Therefore, alternative avenues for
parallelism such as computing multiple timesteps in parallel must be
considered. This has already been attempted almost 60 years ago [5] but only
recently gained greatly increased interest. As a result, a large variety of
so-called parallel-in-time integration algorithms mostly based on Parareal [6]
have been developed in the last decade [7], such as Multigrid Reduction In
Time (MGRIT) [8] or the Parallel Full Approximation Scheme in Space and Time
(PFASST) [9], to name a few. The aforementioned algorithms have been applied
successfully to a variety of problems [10, 11, 12, 13, 14], enhancing parallel
efficiency and allowing codes to scale much further than with spatial
parallelism alone. Nonetheless, parallel-in-time integration has not yet been
adapted to micromagnetics. In this work, we therefore investigate parallel-in-
time integration of the Landau-Lifshitz-Gilbert equation (LLG) with the PFASST
algorithm to gain an understanding of the potential benefits and difficulties.
After explaining the numerical prerequisites and implementation, numerical
experiments giving a first impression on the performance of the scheme are
presented and discussed.

##  II Methods

In the following, basic concepts of micromagnetism and the prerequisites for
the implemented PFASST scheme will be briefly explained, followed by a
description of the algorithm itself. For a more comprehensive overview of
micromagnetics refer to [15] and for more information regarding parallel-in-
time integration, consider visiting the parallel-in-time.org website.

###  II.1 Micromagnetics and the LLG

The micromagnetic model strikes a balance between accuracy and efficiency by
assuming that the quantum mechanical exchange interaction between neighboring
spins is strong enough to align magnetic moments in parallel at distances much
greater than the underlying lattice constant, a characteristic property of
ferromagnets. As a result, the magnetization of a material can be written as a
continuous, vector-valued function
$\bm{M}(\bm{x})=M_{\text{s}}\,\bm{m}(\bm{x})$ where
$\left\lVert\bm{m}(\bm{x})\right\rVert=1$ and $M_{\text{s}}$ refers to the
constant saturation magnetization. Energy contributions influencing the
magnetization are then expressed as integrals involving the magnetization
function. Notable contributions to the total energy $E$ are for example the
Zeeman energy originating from an external field, the demagnetizing energy
resulting from the field caused by the material itself, the aforementioned
exchange energy, or effects such as anisotropy. Furthermore, $E$ is used to
define the so-called effective field $\bm{H}^{\text{eff}}$,

$\bm{H}^{\text{eff}}=-\frac{1}{\mu_{0}\,M_{\text{s}}}\frac{\updelta
E}{\updelta m}\,,$ (1)

where
$\mu_{0}\approx$1.2566\text{\times}{10}^{-6}\text{\,}\mathrm{N}\mathrm{/}\mathrm{A}^{2}$$
refers to the vacuum permeability and $\updelta$ denotes the functional
derivative. While the total energy can be minimized directly to find
equilibrium configurations, the importance of the effective field lies in
modeling magnetization dynamics when inserted into the LLG,

$\frac{\partial\bm{m}}{\partial
t}=-\frac{\gamma}{1+\alpha^{2}}\,\bm{m}\times\bm{H}^{\text{eff}}-\frac{\alpha\,\gamma}{1+\alpha^{2}}\,\bm{m}\times(\bm{m}\times\bm{H}^{\text{eff}})\,,$
(2)

where $\alpha$ refers to a dimensionless damping parameter modeling the
relaxation of the magnetization toward the field lines, whereas the reduced
gyromagnetic ratio
$\gamma=\gamma_{e}/\mu_{0}\approx$2.2128\text{\times}{10}^{5}\text{\,}\mathrm{m}\mathrm{/}{\mathrm{A}\mathrm{s}}$$
determines the frequency of the precession around the field lines. Spatial
discretization of the LLG is commonly done with either finite-difference (FDM)
[1, 2, 3, 4] or finite-element (FEM) [16, 17, 18] methods. In this work, the
finite-element micromagnetic solver magnum.pi [16] is used. It is based on
firedrake [19], which enables distributing meshes and therefore spatial
parallelism for finite-element discretizations of partial differential
equations. By default, time integration in magnum.pi is done with the SUNDIALS
CVODE [20] module. More specifically, a preconditioned, adaptive second-order
Backward Differentiation Formula (BDF) timestepper described in [21] is
employed. This fully implicit time integration method is suited well for
dealing with the high stiffness introduced by the exchange interaction.

###  II.2 Spectral Deferred Corrections

The fundamental ingredient of PFASST is a class of timesteppers referred to as
Spectral Deferred Correction (SDC) methods [22]. Such timesteppers solve
ordinary differential equations by iteratively correcting an initial
approximation using only low order methods. For the sake of brevity, the
method is commonly derived using a preconditioned Richardson iteration
viewpoint [23, 24] as follows: Consider an initial value problem in integral
form,

$u(t)=u_{\text{a}}+\int_{t_{\text{a}}}^{t}f\big{(}\tau,u(\tau)\big{)}\,\mathop{}\\!\mathrm{d}\tau\,,~{}u(t_{\text{a}})=u_{\text{a}}\,,$
(3)

to be integrated until time $t_{\text{b}}$. Constructing an SDC scheme
carrying out this integration begins with the choice of outer quadrature nodes
($t_{1},\dots,t_{m}$) on the interval $[t_{\text{a}},t_{\text{b}}]$ according
to some quadrature rule (e.g. Gauss-Legendre). Then, integrals of the right-
hand side $f$ to these outer quadrature nodes are approximated with the
integration matrix $\bm{Q}\in\mathbb{R}^{m\times m}$:

$\int_{t_{\text{a}}}^{t_{i}}f\big{(}\tau,u(\tau)\big{)}\mathop{}\\!\mathrm{d}\tau\approx\sum_{j=1}^{m}Q_{ij}\,f\big{(}t_{j},u(t_{j})\big{)}\,,\quad
i=1,\dots m$ (4)

In other words, row $i$ of matrix $\bm{Q}$ contains quadrature weights that
integrate from $t_{\text{a}}$ to $t_{i}$. Obtaining said weights usually
involves an interpolation step [25]. Inserting the approximations (4) into (3)
gives rise to the so-called collocation problem

$(\bm{I}_{m}-\bm{Q}\bm{F})[\bm{u}]=\bm{u}_{\text{a}}\,,$ (5)

which, when solved, gives an approximation to the solution at times $t_{i}$,
$\bm{u}=(u_{1},\dots,u_{m})$ where $u_{i}=u(t_{i})$. Here,
$\bm{u}_{\text{a}}=(u_{\text{a}},\dots,u_{\text{a}})\in\mathbb{R}^{m}$ refers
to a vector filled with the initial value, $\bm{I}_{m}$ to the $m$-dimensional
identity operator and
$\bm{F}[\bm{u}]=(f(t_{1},u_{1}),\dots,f(t_{m},u_{m}))\in\mathbb{R}^{m}$
evaluates the right-hand side at the outer quadrature nodes. SDC methods now
solve the nonlinear equation (5) by applying Richardson iteration with a
preconditioner $\bm{Q}^{\text{P}}$:

$(\bm{I}_{m}-\bm{Q}^{\text{P}}\bm{F})[\bm{u}^{k+1}]=(\bm{Q}-\bm{Q}^{\text{P}})\bm{F}[\bm{u}^{k}]+\bm{u}_{\text{a}}$
(6)

The preconditioner $\bm{Q}^{\text{P}}$ usually is a simpler, lower triangular
integration rule for the subintervals, e.g. backward integration. If the
preconditioner is strictly lower triangular (e.g. forward integration),
solving a (non)linear system for the iteration is not required and the scheme
is explicit. Computing one iteration of (6) is often called an SDC sweep.
Carrying out multiple sweeps sequentially improves the accuracy of the
approximation. In fact, one SDC sweep raises the order of accuracy in the
timestep by one or more (depending on the accuracy of the preconditioner), up
until the order of the outer quadrature [26]. Furthermore, using two separate
preconditioners, one of which lower triangular
($\bm{Q}^{\text{P}}_{\text{im}}$) and one strictly lower triangular
($\bm{Q}^{\text{P}}_{\text{ex}}$) allows for constructing an implicit-explicit
splitting (IMEX) version of SDC. With a correspondingly split right-hand side
$\bm{F}=\bm{F}_{\text{ex}}+\bm{F}_{\text{im}}$, the sweeps then look as
follows:

$(\bm{I}_{m}-\bm{Q}^{\text{P}}_{\text{im}}\,\bm{F}_{\text{im}}-\bm{Q}^{\text{P}}_{\text{ex}}\,\bm{F}_{\text{ex}})[\bm{u}^{k+1}]=\\\
(\bm{Q}\bm{F}-\bm{Q}^{\text{P}}_{\text{im}}\,\bm{F}_{\text{im}}-\bm{Q}^{\text{P}}_{\text{ex}}\,\bm{F}_{\text{ex}})[\bm{u}^{k}]+\bm{u}_{\text{a}}$
(7)

###  II.3 Multilevel SDC

Multigrid methods outsource computationally expensive tasks to a coarser
discretization where they can be carried out in a cheap manner. This can be
used to accelerate convergence for SDC timesteppers and is used in PFASST to
facilitate a cheap estimation of the solution at future timesteps. In so-
called Multilevel SDC (MLSDC) methods, the Full Approximation Scheme (FAS)
multigrid correction technique [27] is used to obtain a correction for the
right-hand side integrals at the coarse level,

$\bm{\tau}=\bm{T}_{\text{f}}^{\text{c}}\,\bm{Q}\,\bm{F}[\bm{u}^{k}]-\tilde{\bm{Q}}\,\bm{T}_{\text{f}}^{\text{c}}\,\bm{F}[\bm{u}^{k}]\,,$
(8)

where $\bm{T}_{\text{f}}^{\text{c}}$ refers to a restriction operator which
can transform onto both a coarser mesh and time discretization (i.e. fewer
nodes in the outer quadrature). Furthermore, $\tilde{\bm{Q}}$ represents the
integration matrix on the coarse level. Adding $\bm{\tau}$ to the right-hand
side of (6) for the coarse level SDC sweeps enhances accuracy, often described
as computing the solution at the resolution of the coarse grid but with the
accuracy of the fine grid. Overall, a complete MLSDC iteration involves the
following:

  1. 1.

Restrict the fine level approximation:
$\tilde{\bm{u}}^{k}=\bm{T}_{\text{f}}^{\text{c}}[\bm{u}^{k}]$.

  2. 2.

Compute $\bm{\tau}$ according to (8) and carry out corrected SDC sweeps on the
coarse level to obtain $\tilde{\bm{u}}^{k+1/2}$.

  3. 3.

Get a coarse level correction for the solution:
$\tilde{\bm{\delta}}^{k}=\tilde{\bm{u}}^{k+1/2}-\tilde{\bm{u}}^{k}$.

  4. 4.

Prolong the correction to the fine level with the prolongation operator
$\bm{T}_{\text{c}}^{\text{f}}$,
$\bm{\delta}^{k}=\bm{T}_{\text{c}}^{\text{f}}[\tilde{\bm{\delta}}^{k}]$, and
apply it to the fine level approximation:
$\bm{u}^{k+1/2}=\bm{u}^{k}+\bm{\delta}^{k}$.

  5. 5.

Carry out fine level SDC sweeps to obtain $\bm{u}^{k+1}$.

Aside from the two-level example above, this can also be carried out on three
or more levels [9, 25].

###  II.4 Parallel Full Approximation Scheme in Space and Time

The structure of the MLSDC algorithm shown above is suited well for
introducing parallelism in time. All that remains to do is strategically
passing intermediate results forward, to be used as initial values for
processes that calculate future timestep values. Consider Fig. 1, which
depicts two-level MLSDC iterations that take place on multiple processes
($P_{0},\dots,P_{3}$), computing multiple timesteps
($[t_{0},t_{1}],[t_{1},t_{2}],\dots$) concurrently. Aside from an
initialization procedure at the beginning, the only additions to the MLSDC
algorithm are passing the coarse level result forward to the next process upon
completing coarse sweeps and passing fine level results forward to the next
iteration on the next process. Regarding the initialization procedure, the
repeated coarse sweeps shown in the Fig. are a common choice.

Figure 1: Two-level PFASST iterations on four processes visualized with
pfasst-tikz [28]. End results of coarse (blue) and fine (red) sweeps are sent
forward to the next process (arrows), where they are used as initial values.

There are many variants of PFASST with different communication patterns or
other modifications. The pySDC library [29] for example, which we used to
prototype the PFASST and SDC solvers implemented in this work, adopts a
multigrid viewpoint. While this enables a more general implementation of the
algorithm, it also decreases efficiency. Here however we use the classical
two-level variant depicted in Fig. 1. When compared to serial SDC, the
parallel efficiency of this variant is at most $K_{\text{s}}/K_{\text{p}}$
[30], where $K_{\text{s}}$ refers to the number of iterations required to
solve the problem with SDC, whereas $K_{\text{p}}$ refers to the required
PFASST iterations to solve the problem with the same accuracy. This is an
improvement over the upper efficiency bound of $1/K$ for Parareal [6], where
$K$ is the number of iterations required for convergence.

###  II.5 Implementation

The first of many choices in implementing a PFASST scheme is choosing a
numerical quadrature for the SDC sweeps. Usually, quadrature types with the
highest possible order given a certain number of nodes are chosen. Gauss-
Legendre quadrature for instance integrates polynomials of up to degree
$2\,n-1$ exactly with $n$ nodes and is optimal in this sense. On the other
hand, Gauss-Legendre nodes on a given interval do not include the endpoint of
said interval, where the desired solution might lie. Hence, an additional
routine for calculating the endpoint is needed to integrate ODEs with such a
quadrature type. This issue can be avoided entirely with quadratures that
include the endpoint, such as Gauss-Radau quadrature which integrates
polynomials of up to degree $2\,n-2$ exactly with $n$ nodes. In this work,
both Gauss-Legendre and Gauss-Radau nodes are used, where the solution at the
endpoint is computed by integrating with the outer quadrature weights if
necessary.

Next, the preconditioners must be assembled. Common choices here are forward
integration as an explicit preconditioner and the so-called LU-trick [31],

$\bm{Q}^{\text{P}}_{\text{LU}}=\bm{U}^{\text{T}}\,,\text{ where
}\bm{Q}^{\text{T}}=\bm{L}\,\bm{U}\,,$ (9)

as an implicit preconditioner which performs favorably for stiff problems.
Both preconditioning methods are used in the code, because the LLG is not only
stiff, but also a good candidate for implicit-explicit time integration: On
the one hand, the stiffness is mainly introduced by the short-ranged exchange
interaction, which is relatively cheap to compute and can thus be integrated
implicitly without issue. On the other hand, the LLG permits a very
straightforward splitting of the right-hand side into the contributions of
individual field terms, as it is linear in the effective field. In the
implemented code, users can choose which field terms are integrated explicitly
or implicitly. One major benefit of using such an IMEX time integrator for the
LLG is that repeated evaluations of the expensive demagnetizing field can be
avoided. Therefore, we always use the IMEX variant of SDC in the numerical
experiments of section III, where only the exchange field is evaluated
implicitly.

The nonlinear system that arises in the SDC sweeps is solved with line search
Newton iteration, more specifically the PETSc [32] implementation with the
Flexible Generalized Minimum Residual Method (FGMRES) as the linear solver.
The necessary right-hand side and Jacobian evaluation routines are already
part of magnum.pi. Together with firedrake, this enables spatial parallelism.

Finally, transfer operators between the coarse and fine space/time
discretizations must be implemented. In this work, the firedrake.multigrid
module is used for spatial transfers, as it allows for restriction and
prolongation between distributed meshes. For the temporal domain, restriction
and prolongation are realized with interpolation between the coarse and fine
level quadrature nodes.

Regarding error control, serial SDC timesteppers have access to a very simple
error estimate: Since the order of accuracy of the solution increases with
every sweep, subtracting the solution values of the previous and current
iteration gives an estimate for the local truncation error of the previous
iterate. The implemented stepsize adaptation is then straightforward: If the
error tolerance is not met, the stepsize is decreased by $50\%$, whereas the
stepsize is increased by $25\%$ after 10 accepted steps. Additionally, as long
as the error exceeds $80\%$ of the tolerance, the stepsize is not increased
further. As an alternative to the error estimate, the maximum residual of the
collocation problem (5) can be considered for error control. This is
particularly useful for comparing the accuracy of SDC and PFASST results.

Lastly, it is worth mentioning that the implemented PFASST controller is very
similar in structure to the classical two-level variant in Fig. 1, with some
modifications for increased efficiency. For instance, reevaluating function
values when returning from the coarse to the fine grid can be avoided by also
correcting the function values in the same way as the approximate solution.
Also, CPU idle times are avoided as much as possible by immediately starting
computation on a new timestep once a timestep has converged, instead of
waiting for a complete set of timesteps to converge and restarting the PFASST
algorithm to compute the next set.

##  III Numerical Experiments

In the following, we demonstrate the efficiency of the implemented
timestepping methods by simulating various micromagnetic processes. Exploring
not only novel parallel timestepping methods but also new kinds of High
Performance Computing (HPC) infrastructure, simulations are carried out on the
Google Cloud. Here, the HPC toolkit [33] is used to first build a magnum.pi
Virtual Machine (VM) image and then create single VM instances or even entire
SLURM [34] clusters on the cloud. In this work, the c2d machine type is used
for all experiments. More specifically, unless stated otherwise simulations
are carried out on c2d-highcpu-32 nodes with 16 physical processor cores each,
where the compact placement rule for better inter-node communication has been
enabled. For another example of micromagnetics code running on cloud services,
consider [35].

###  III.1 Standard Problem 4

PFASST relies entirely on SDC methods for timestepping. As such, investigating
the performance of serial SDC timestepping methods gives valuable insight
about the potential performance of a PFASST scheme. This preliminary
experiment examines SDC as a serial timestepping method for the LLG by
comparing it to the BDF method normally employed in magnum.pi. Additionally, a
first small-scale time-parallel run is attempted. As parallelism in space is
not of interest here, we simulate the well-known $\mu$MAG Standard Problem #4
(SP4) [36].

First, we simulate SP4 with magnum.pi’s built-in adaptive BDF solver and the
adaptive SDC timestepper at varying tolerances. These runs are compared to a
reference solution obtained with BDF at very strict ($tol$ = ${10}^{-10}$)
error tolerances (Tab. 1, Fig. 2).

Table 1: Tolerances $tol$, runtimes $t$, average relative $L^{2}$-norm errors $\overline{e}_{\text{rel}}$ and number of external $n_{\text{ext}}$, demagnetizing $n_{\text{dem}}$ as well as exchange $n_{\text{exc}}$ field evaluations for SP4. | $tol$ |  $t$ (s) | $\overline{e}_{\text{rel}}$ | $n_{\text{ext}}$ | $n_{\text{dem}}$ | $n_{\text{exc}}$
---|---|---|---|---|---|---
SDC | $1.0\text{\times}{10}^{-2}$ | $47$ | $1.54\text{\times}{10}^{-3}$ | 2324 | 6426
BDF | $1.0\text{\times}{10}^{-5}$ | $66$ | $3.42\text{\times}{10}^{-3}$ | 5124
SDC | $1.0\text{\times}{10}^{-5}$ | $140$ | $9.66\text{\times}{10}^{-6}$ | 6111 | 27798
BDF | $2.5\text{\times}{10}^{-9}$ | $708$ | $1.20\text{\times}{10}^{-5}$ | 67587
Figure 2: Average magnetization in the coordinate directions for SP4. The BDF
run with $tol=${10}^{-10}$$, SDC run with $tol=${10}^{-5}$$ and PFASST run
with $n_{\text{t}}=3$ are shown.

On average, the SDC runs are slightly more accurate than comparable BDF runs
but also significantly faster, especially at high accuracy. The speedup is
explained as follows: First, a large fraction of the computationally costly
stray field evaluations is avoided by employing the IMEX variant of SDC
instead of the fully implicit BDF method. Additionally, two Gauss-Legendre
nodes were used for the high-accuracy SDC run. This results in a third order
timestepping method as opposed to the second order BDF method.

Next, we determine the error of individual SDC iterations on the first
timestep of SP4. To this end, we use five Gauss-Legendre nodes for SDC and a
BDF reference solution with even stricter tolerances ($tol=${10}^{-14}$$).
Plotting the errors obtained this way against the iteration number confirms
the theoretically predicted exponential convergence (Fig. 3). Also, one can
clearly see how convergence stops after the accuracy of the quadrature is
exhausted and how the error estimate gives a good approximation of the error
until then.

Figure 3: Relative $L^{2}$-norm error $e_{\text{rel}}$ and error estimate at
SDC iteration $n$ on the first timestep of SP4.

Finally, since the PFASST controller does not yet support adaptive stepsizes,
we use a fixed stepsize of $dt=${10}^{-12}\text{\,}\mathrm{s}$$ and a relative
residual tolerance of $restol=${10}^{-6}$$ to compute SP4 parallel-in-time on
up to three cores (Tab. 2). Here, the fine level mesh is the same as in the
previous experiments with 15000 elements in total, whereas the coarse level
mesh consists of only 600 elements. The residual tolerance is chosen such that
the results roughly match the accuracy of the SDC/BDF runs in Tab. 1.

Table 2: Number of cores $n_{\text{t}}$, runtimes $t$, speedups $S$ and average relative $L^{2}$-norm errors $\overline{e}_{\text{rel}}$ when computing SP4 with PFASST. | $n_{\text{t}}$ |  $t$ (s) | $S$ | $\overline{e}_{\text{rel}}$
---|---|---|---|---
PFASST | 1 | $246$ | 1.00 | $3.87\text{\times}{10}^{-6}$
2 | $148$ | 1.66 | $5.29\text{\times}{10}^{-6}$
3 | $114$ | 2.16 | $1.27\text{\times}{10}^{-5}$

As can be seen in Tab. 2, using PFASST on a single core (i.e. MLSDC) is
initially slower than using adaptive SDC. This relatively large slowdown is
primarily caused by the already very small fine mesh size, which makes the
additional overhead introduced in PFASST much more noticeable and also
diminishes the speedup gained from transferring to the coarse level. Despite
this, a parallel speedup of up to 2.16 is achieved with three cores in time.
Note that the average relative error increases with the number of cores in
time. However, this is expected and originates from the fact that one single
PFASST iteration can significantly reduce the relative residual. Therefore,
while all results meet the same residual tolerance, they vary in accuracy
depending on how much the relative residual decreases below the tolerance on
the final iteration of every timestep. It is also worth mentioning that time
parallelism on more than a few cores is likely not very efficient for problems
like SP4, unless very high accuracies are desired. This is due to the fact
that the magnetization changes quickly in SP4, making future timesteps hard to
predict. Consequently, parallel efficiency decreases disproportionately as
more and more cores are used and the time interval that is computed in
parallel grows. Nonetheless, for the specified residual tolerance the time-
parallel PFASST runs already offer a small speedup over serial runs.

###  III.2 Airbox Hysteresis

Unlike the switching process in the previous subsection, magnetization
generally changes slowly when computing hysteresis loops. This should improve
parallel efficiency and allow using more cores for time parallelism.
Therefore, we investigate this problem type next: We set a
$$1000\text{\,}\mathrm{nm}$\times$1000\text{\,}\mathrm{nm}$\times$20\text{\,}\mathrm{nm}$$
rectangular domain to a constant initial magnetization of
$\bm{m}=(0.01,0,-1)^{\text{T}}$ before normalization and then reverse the
magnetization by simulating the LLG for
$t_{\text{end}}=$2\text{\,}\mathrm{ns}$$ with an external field of the form

$\bm{H}_{\text{ext}}(t)=\left(\frac{2t}{t_{\text{end}}}-1\right)\,\frac{B_{\text{max}}}{\mu_{0}}\,\bm{e}_{z}\,,~{}~{}B_{\text{max}}=$1\text{\,}\mathrm{T}$\,,$

where $\bm{e}_{\text{z}}$ refers to the unit vector in $z$-direction. The
material parameters are the same as in SP4, except for the damping where we
set $\alpha=1$. In other words, we compute one branch of a hysteresis loop in
the direction of the hard axis. Because the hybrid finite-element-boundary-
element (FEM-BEM) demagnetizing field solver in magnum.pi does not yet support
spatial parallelism, this is done with the airbox solver module and a
surrounding airbox of the dimensions
$$4000\text{\,}\mathrm{nm}$\times$4000\text{\,}\mathrm{nm}$\times$400\text{\,}\mathrm{nm}$$.
We use a mesh with 33033 elements on the coarse PFASST level and, after
refinement by bisection, 264264 elements on the fine level. With seven Gauss-
Radau nodes on the fine and coarse level, a fixed stepsize of
$dt=$2\text{\times}{10}^{-11}\text{\,}\mathrm{s}$$ and a relative residual
tolerance of $restol=${10}^{-8}$$, we obtain the runtimes in Fig. 4 on a
single node. To ensure correctness, we also carry out a BDF run with
$tol=${10}^{-7}$$ and plot the hysteresis curve against selected PFASST
results (Fig. 5).

Figure 4: PFASST timings for the hysteresis problem, solved with
$n_{\text{s}}$ cores in space and $n_{\text{t}}$ parallel timesteps on a
single node. Figure 5: Hysteresis curves computed with BDF and PFASST, where
$n_{\text{t}}$ time intervals are computed in parallel and $H_{\text{ext}}$
refers to the $z$-component of the external field.

Notably, up until eight cores, pure time parallelism is more efficient than
pure spatial parallelism. This is most likely a consequence of the still
relatively small mesh used in the simulations. On the other hand, the
efficiency of the time-parallel approach drops quite quickly, such that using
16 cores in time actually slows the computation down compared to using only
eight cores.

Finally, Fig. 6 compares the scaling of space-parallel only runs with space-
and time-parallel runs on one (solid lines) or several (dashed lines) nodes.

Figure 6: Strong scaling runtimes $t$ for the hysteresis problem, utilizing
$n_{\text{tot}}$ cores on individual nodes (solid lines) and across several
nodes (dashed lines). The ideal scaling (dashed-dotted line) is inversely
proportional to $n_{\text{tot}}$.

Overall, by using 128 cores in total the runtime could be reduced by a factor
of 20.56, where the speedup of 5.75 with spatial parallelism alone is enhanced
by an additional factor of 3.58 through time parallelism. One should however
note that the residual tolerance used in this numerical experiment is
relatively strict for simulating such a simple problem, which increases the
parallel efficiency of PFASST. On top of that, the coarse level resolution is
already quite good, which similarly benefits PFASST. In any case, the results
show that, unlike spatial parallelism, parallel-in-time integration can be
used efficiently across multiple nodes. This observation is especially
interesting in the context of micromagnetic GPU codes, where time parallelism
could offer a straightforward way to implement multi-GPU extensions.

###  III.3 Magnonic damping

Finally, we consider a problem of oscillatory nature, where high accuracies
are often desired. In [37], a spin pumping model for magnum.pi based on spin-
diffusion is implemented and used to simulate the propagation of a magnon
damped by spin torque. More specifically, the magnon is driven by an
oscillating magnetic field on one side of a
$$1251\text{\,}\mathrm{nm}$\times$10\text{\,}\mathrm{nm}$\times$10\text{\,}\mathrm{nm}$$
waveguide with uniaxial anisotropy along the $y$-axis and damped by a spin
sink placed on the waveguide starting from $y=$251\text{\,}\mathrm{nm}$$.
Furthermore, the simulation is carried out with zero damping ($\alpha=0$) and
with the demagnetizing field disabled. We simulate this magnonic damping
process for $3\text{\,}\mathrm{ns}$ with the same waveguide structure and
material parameters. As for the solver parameters, the PFASST runs utilize a
coarse mesh with 84436 elements and a fine mesh resulting from bisection with
675488 elements alongside three Gauss-Radau nodes, a fixed stepsize of
$dt=$5\text{\times}{10}^{-12}\text{\,}\mathrm{s}$$ as well as a relative
residual tolerance of $rtol=${10}^{-5}$$. Because the spin torque calculations
are relatively expensive on a mesh of this size, we now use the
c2d-highcpu-112 machine type on the Google Cloud, with 56 physical cores per
node. Once again, we first simulate on a single node and then extend the
parallelism by computing timesteps in parallel on up to four nodes, for a
total maximum of 224 cores. This results in the strong scaling depicted in
Fig. 7.

Figure 7: Strong scaling runtimes $t$ for the magnon damping problem,
utilizing $n_{\text{tot}}$ cores on individual nodes (solid lines) and across
several nodes (dashed lines). The ideal scaling (dashed-dotted line) is
inversely proportional to $n_{\text{tot}}$.

At 56 cores, a speedup of 17.56 is achieved with pure spatial parallelism and
an additional speedup of 2.08 when using four such nodes parallel-in-time for
a total speedup of 36.51. As comparing the performance of PFASST to BDF is
also interesting here, we compute a BDF solution with similar accuracy. For
this experiment, the quantity of interest is the magnitude of the transversal
magnetization component along the $y$-axis, $m_{\perp}$, at
$t=$2.4\text{\,}\mathrm{ns}$$ (Fig. 8). Therefore, we use this quantity to
estimate the accuracies of our solutions by comparing to a reference BDF run
with $tol=${10}^{-7}$$ (Tab. 3).

Figure 8: Transversal magnetization component magnitude along the $y$ axis $m_{\perp}$ at $t=$2.4\text{\,}\mathrm{ns}$$ in the magnon damping simulations for selected time integration methods and tolerances. Table 3: Number of cores used in space $n_{\text{s}}$, parallel timesteps $n_{\text{t}}$, runtimes $t$ and average relative errors $\overline{e}_{\text{rel}}$ of $m_{\perp}$ for the magnon damping problem. | $n_{\text{s}}$ | $n_{\text{t}}$ |  $t$ (h) | $\overline{e}_{\text{rel}}$
---|---|---|---|---
BDF ($tol=${10}^{-5}$$) | 1 | 1 | 36.11 | $6.12\text{\times}{10}^{-04}$
PFASST | 1 | 1 | 17.14 | $4.24\text{\times}{10}^{-04}$
14 | 1 | 2.63 | $4.24\text{\times}{10}^{-04}$
28 | 1 | 1.70 | $5.00\text{\times}{10}^{-04}$
14 | 2 | 1.60 | $4.75\text{\times}{10}^{-04}$
56 | 1 | 0.98 | $4.23\text{\times}{10}^{-04}$
14 | 4 | 0.94 | $7.61\text{\times}{10}^{-04}$
28 | 2 | 0.81 | $5.40\text{\times}{10}^{-04}$
56 | 2 | 0.67 | $4.74\text{\times}{10}^{-04}$
56 | 3 | 0.55 | $3.16\text{\times}{10}^{-04}$
56 | 4 | 0.47 | $7.60\text{\times}{10}^{-04}$

As explained in subsection III.1, although all PFASST runs meet the same
residual tolerance, the error generally increases when integrating parallel-
in-time. It is also worth noting again here that all field terms except for
the exchange field are integrated explicitly when using PFASST. This includes
the expensive spin torque evaluations and already roughly halves the runtime
compared to BDF on a single core. Lastly, as opposed to the hysteresis
simulation in the previous chapter, the error tolerances chosen for this
numerical experiment are on the lower end. Higher accuracies are generally
desired in this context, where PFASST should perform even better overall.

##  IV Conclusions and Discussion

In this work, we implemented IMEX SDC/PFASST timesteppers for the LLG as
extensions to the micromagnetic FEM solver magnum.pi. First numerical
experiments lead to the following observations:

IMEX Spectral Deferred Correction method timesteppers are competitive with
fully implicit preconditioned BDF(2) methods at moderate accuracies. At higher
accuracies, the flexibility of SDC methods enables straightforward
construction of high-order IMEX timestepping methods for greatly improved
performance.

Parallel-in-time integration of the LLG with decent efficiencies has been
achieved across a variety of problem classes (Tab. 4), where the code
generally performs better when high accuracies are desired. Combining space-
and time-parallelism leads to better efficiencies than using only space-
parallelism — in particular, parallel-in-time integration is also suited well
for parallelism across multiple nodes, where it was successfully used to scale
beyond spatial limits.

Table 4: Maximum speedups $S$ over serial PFASST in the numerical experiments and corresponding parallel efficiencies $\eta=S/n_{\text{tot}}$, where $n_{\text{tot}}$ refers to the total number of cores used. | SP4 | Hysteresis | Magnon
---|---|---|---
| $S$ | $\eta$ | $S$ | $\eta$ | $S$ | $\eta$
Space | - | - | 5.75 | 0.36 | 17.56 | 0.31
Time | 2.16 | 0.72 | 3.58 | 0.45 | 2.08 | 0.52
Total | 2.16 | 0.72 | 20.56 | 0.16 | 36.51 | 0.16

These preliminary results are best understood as a proof of concept for
parallel-in-time integration in micromagnetic simulations, since many
optimizations and improvements are yet to be made. For instance, further
optimizing the linear and nonlinear solvers leads to faster convergence and
enables larger timesteps, which would especially increase performance at low
to moderate accuracies. Furthermore, enhancements to SDC such as linearly
implicit formulations or variants like inexact SDC [38] or Krylov-SDC [24] can
yield better convergence. Options for additionally parallelizing SDC “across
the method” are also available [39] and the adaptiveness of the algorithm
could be improved by increasing/decreasing the number of quadrature nodes
automatically as needed. Similarly, adaptive timesteps could be implemented
for the PFASST controller. There are also many other parallel-in-time
integration algorithms besides PFASST that have yet to be tested in the
context of integrating the LLG. Finally, extending GPU simulation codes to
multiple GPUs with time parallelism seems promising due to the relatively
straightforward implementation. This and the above improvements will be
investigated in the future.

Acknowledgements

This research was funded in whole, or in part, by the Austrian Science Fund
(FWF) P 34671. For the purpose of open access, the author has applied a CC BY
public copyright license to any Author Accepted Manuscript version arising
from this submission.

## References

  * Donahue [1999] M. Donahue,  _OOMMF User’s Guide, Version 1.0_ (1999).
  * Vansteenkiste  _et al._ [2014] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, The design and verification of MuMax3, AIP Advances 4, 107133 (2014).
  * Bruckner  _et al._ [2023] F. Bruckner, S. Koraltan, C. Abert, and D. Suess, magnum.np: a PyTorch based GPU enhanced finite difference micromagnetic simulation framework for high level development and inverse design, Scientific Reports 13, 12054 (2023).
  * Lepadatu [2023] S. Lepadatu,  _Accelerating micromagnetic and atomistic simulations using multiple GPUs_, Tech. Rep. (2023) arXiv:2308.08447 [cond-mat, physics:physics] type: article.
  * Nievergelt [1964] J. Nievergelt, Parallel methods for integrating ordinary differential equations, Communications of the ACM 7, 731 (1964).
  * Lions  _et al._ [2001] J.-L. Lions, Y. Maday, and G. Turinici, Résolution d’EDP par un schéma en temps «pararéel », Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 332, 661 (2001).
  * Gander [2015] M. J. Gander, 50 Years of Time Parallel Time Integration, in  _Multiple Shooting and Time Domain Decomposition Methods_, Contributions in Mathematical and Computational Sciences, edited by T. Carraro, M. Geiger, S. Körkel, and R. Rannacher (Springer International Publishing, Cham, 2015) pp. 69–113.
  * Friedhoff  _et al._ [2012] S. Friedhoff, R. D. Falgout, T. V. Kolev, S. MacLachlan, and J. B. Schroder,  _A Multigrid-in-Time Algorithm for Solving Evolution Equations in Parallel_, Tech. Rep. (2012).
  * Emmett and Minion [2012] M. Emmett and M. Minion, Toward an efficient parallel in time method for partial differential equations, Communications in Applied Mathematics and Computational Science 7, 105 (2012).
  * Fischer  _et al._ [2005] P. F. Fischer, F. Hecht, and Y. Maday, A Parareal in Time Semi-implicit Approximation of the Navier-Stokes Equations, in  _Domain Decomposition Methods in Science and Engineering_, Lecture Notes in Computational Science and Engineering, edited by T. J. Barth, M. Griebel, D. E. Keyes, R. M. Nieminen, D. Roose, T. Schlick, R. Kornhuber, R. Hoppe, J. Périaux, O. Pironneau, O. Widlund, and J. Xu (Springer, Berlin, Heidelberg, 2005) pp. 433–440.
  * Legoll  _et al._ [2022] F. Legoll, T. Lelièvre, and U. Sharma, An adaptive parareal algorithm: Application to the simulation of molecular dynamics trajectories, SIAM Journal on Scientific Computing 44, B146 (2022), https://doi.org/10.1137/21M1412979 .
  * Schöps  _et al._ [2018] S. Schöps, I. Niyonzima, and M. Clemens, Parallel-In-Time Simulation of Eddy Current Problems Using Parareal, IEEE Transactions on Magnetics 54, 1 (2018).
  * Speck  _et al._ [2014] R. Speck, D. Ruprecht, R. Krause, M. Emmett, M. Minion, M. Winkel, and P. Gibbon, Integrating an N-Body Problem with SDC and PFASST, in  _Domain Decomposition Methods in Science and Engineering XXI_, Lecture Notes in Computational Science and Engineering, edited by J. Erhel, M. J. Gander, L. Halpern, G. Pichot, T. Sassi, and O. Widlund (Springer International Publishing, Cham, 2014) pp. 637–645.
  * Götschel and Minion [2018] S. Götschel and M. L. Minion, Parallel-in-Time for Parabolic Optimal Control Problems Using PFASST, in  _Domain Decomposition Methods in Science and Engineering XXIV_, Lecture Notes in Computational Science and Engineering, edited by P. E. Bjørstad, S. C. Brenner, L. Halpern, H. H. Kim, R. Kornhuber, T. Rahman, and O. B. Widlund (Springer International Publishing, Cham, 2018) pp. 363–371.
  * Abert [2019] C. Abert, Micromagnetics and spintronics: models and numerical methods, The European Physical Journal B 92, 120 (2019).
  * Abert  _et al._ [2013] C. Abert, L. Exl, F. Bruckner, A. Drews, and D. Suess, magnum.fe: A micromagnetic finite-element simulation code based on FEniCS, Journal of Magnetism and Magnetic Materials 345, 29 (2013).
  * Chang  _et al._ [2011] R. Chang, S. Li, M. V. Lubarda, B. Livshitz, and V. Lomakin, FastMag: Fast micromagnetic simulator for complex magnetic structures (invited), Journal of Applied Physics 109, 07D358 (2011).
  * Körber  _et al._ [2022] L. Körber, G. Quasebarth, A. Hempel, F. Zahn, A. Otto, E. Westphal, R. Hertel, and A. Kakay, TetraX: Finite-Element Micromagnetic-Modeling Package (2022).
  * Rathgeber  _et al._ [2016] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. Mcrae, G.-T. Bercea, G. R. Markall, and P. H. J. Kelly, Firedrake: Automating the Finite Element Method by Composing Abstractions, ACM Transactions on Mathematical Software 43, 24:1 (2016).
  * Hindmarsh  _et al._ [2005] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward, SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers, ACM Transactions on Mathematical Software 31, 363 (2005).
  * Suess  _et al._ [2002] D. Suess, V. Tsiantos, T. Schrefl, J. Fidler, W. Scholz, H. Forster, R. Dittrich, and J. J. Miles, Time resolved micromagnetics using a preconditioned time integration method, Journal of Magnetism and Magnetic Materials 248, 298 (2002).
  * Dutt  _et al._ [2000] A. Dutt, L. Greengard, and V. Rokhlin, Spectral Deferred Correction Methods for Ordinary Differential Equations, BIT Numerical Mathematics 40, 241 (2000).
  * Ruprecht and Speck [2016] D. Ruprecht and R. Speck, Spectral deferred corrections with fast-wave slow-wave splitting, SIAM Journal on Scientific Computing 38, A2535 (2016), https://doi.org/10.1137/16M1060078 .
  * Huang  _et al._ [2006] J. Huang, J. Jia, and M. Minion, Accelerating the convergence of spectral deferred correction methods, Journal of Computational Physics 214, 633 (2006).
  * Bolten  _et al._ [2017] M. Bolten, D. Moser, and R. Speck, A multigrid perspective on the parallel full approximation scheme in space and time, Numerical Linear Algebra with Applications 24, e2110 (2017).
  * Causley and Seal [2019] M. F. Causley and D. C. Seal, On the convergence of spectral deferred correction methods, Communications in Applied Mathematics and Computational Science 14, 33 (2019), arXiv:1706.06245 [math].
  * Henson [2002] V. E. Henson,  _Multigrid Methods for Nonlinear Problems: An Overview_, Tech. Rep. (2002).
  * Koehler [2015] F. Koehler, Pfasst tikz, https://github.com/f-koehler/pfasst-tikz/tree/master (2015).
  * Speck [2019] R. Speck, Algorithm 997: pySDC—Prototyping Spectral Deferred Corrections, ACM Transactions on Mathematical Software 45, 35:1 (2019).
  * Minion [2011] M. Minion, A hybrid parareal spectral deferred corrections method, Communications in Applied Mathematics and Computational Science 5, 265 (2011).
  * Weiser [2015] M. Weiser, Faster SDC convergence on non-equidistant grids by DIRK sweeps, BIT Numerical Mathematics 55, 1219 (2015).
  * Balay  _et al._ [2019] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, _et al._ ,  _PETSc users manual_ (2019).
  * Google Cloud Platform [2023] Google Cloud Platform, Google HPC-Toolkit, https://github.com/GoogleCloudPlatform/hpc-toolkit (2023).
  * Yoo  _et al._ [2003] A. B. Yoo, M. A. Jette, and M. Grondona, SLURM: Simple Linux Utility for Resource Management, in  _Job Scheduling Strategies for Parallel Processing_, Lecture Notes in Computer Science, edited by D. Feitelson, L. Rudolph, and U. Schwiegelshohn (Springer, Berlin, Heidelberg, 2003) pp. 44–60.
  * Jermain  _et al._ [2016] C. L. Jermain, G. E. Rowlands, R. A. Buhrman, and D. C. Ralph, GPU-accelerated micromagnetic simulations using cloud computing, Journal of Magnetism and Magnetic Materials 401, 320 (2016).
  * Eicke [2000] J. Eicke, µMag Standard Problem #4 specification, https://www.ctcms.nist.gov/~rdm/std4/spec4.html (2000).
  * Gattringer  _et al._ [2022] M. Gattringer, C. Abert, F. Bruckner, A. Chumak, and D. Suess, Micromagnetically integrated numerical model of spin pumping based on spin diffusion, Physical Review B 106, 10.1103/PhysRevB.106.024417 (2022).
  * Speck  _et al._ [2016] R. Speck, D. Ruprecht, M. Minion, M. Emmett, and R. Krause, Inexact Spectral Deferred Corrections, in  _Domain Decomposition Methods in Science and Engineering XXII_, Lecture Notes in Computational Science and Engineering, edited by T. Dickopf, M. J. Gander, L. Halpern, R. Krause, and L. F. Pavarino (Springer International Publishing, Cham, 2016) pp. 389–396.
  * Speck [2018] R. Speck, Parallelizing spectral deferred corrections across the method, Computing and Visualization in Science 19, 75 (2018).