Sumsets and fixed points of substitutions
F. Michel Dekking

DIAM,  Delft University of Technology, Faculty EEMCS, Delft, The Netherlands.

In this paper we introduce a technique to determine the sumset $A+A$, where $A$ is the indicator function of the 0's occurring in a fixed point $x$ of a substitution on the alphabet $\{0,1\}$.

§ INTRODUCTION

Let $A$ be the sequence given by $A(n)= \lfloor n\alpha\rfloor$ for $n\ge 1$, where $\alpha = (1+\sqrt{5})/2$ is the golden mean. In the recent paper [7], it was proved that the sumset $A+A$ is equal to the natural numbers with exception of the numbers 1 and 3.
Here the sumset $A+A$ is defined by

$$A+A = \{a+b: a\in A, b\in A\}.$$

It was also shown that if $B$ is defined by $B(n)= \lfloor n\alpha^2\rfloor$, then  the set $B+B$ has a complement that is infinite. The determination of $B+B$ takes more than 8 pages, and the  complicated structure of $B+B$ is described as containing “some fractal patterns". The goal of the present paper is to elucidate what the source of these fractal patterns is, and at the same time to introduce a vast generalization of these types of sumsets.

Interestingly, Shallit in [9] also positioned the sumset problem of the paper [7] in the combinatorics on words context, but in a completely different way. The focus there is on the Fibonacci representation (also known as Zeckendorf representation) of the natural numbers, and the proofs are computer assisted. See also the paper [11].

The crucial observation that we make is that the sequences $A=(\lfloor n\alpha\rfloor) =(1, 3, 4, 6, 8, 9, 11,\dots)$ and $B=(\lfloor n\alpha^2\rfloor)=(2, 5, 7, 10, 13, 15,\dots)$  give the positions of 0's, respectively the positions of 1's in the infinite Fibonacci word $x_{\rm F}=010010100100\dots$, fixed point of the substitution $0\mapsto 01, \, 1\mapsto 0$ (see, e.g., [8]). Let $\sigma$ be any substitution on the monoid $\{0,1\}^*$, admitting a fixed point $x=\sigma(x)$. Then, let
$A$ give the positions of 0 in $x$, and let $B$ give the positions of 1 in $x$.

Problem: determine $A+A$, $A+B$ and $B+B$.

The proofs in [7] are purely arithmetical. Our approach, inspired by [4], is to analyse  $A+A$ by studying the product set $A\times A$.
This amounts to passing from fixed points  of substitutions to fixed points of two-dimensional substitutions.  See, e.g., [6] for an overview of the theory of two-dimensional substitutions. Our application of two-dimensional substitutions to the sumset problem is self-contained.

In Section <ref> we give a very simple proof of the  $A+A$  result of [7] using a
two-dimensional substitution which has $x_{\rm F}\times x_{\rm F}$ as fixed point.
A proof of the $B+B$ result from [7], analogously to the way this is done in  Theorem <ref>, is very complex. The numbers of columns on the right side of  the 2D fixed point of the product substitution that one has to use is unbounded for the Fibonacci sunstitution, whereas it is 15 for the von Neumann substitution.

In Section <ref> we solve the $A+A$, $A+B$ and $B+B$ problem for the Thue-Morse  word $t=0110100110010110\dots$, fixed point of  the substitution $0\mapsto 01, \, 1\mapsto 10$.
We remark that our result is equivalent to Theorem 2 of the paper [11], which is proved in a completely different, computer assisted way.

In Section <ref> we solve the $A+A$ and $B+B$ problem for the von Neumann word.
The interest of this example lies in the fact that this example is one of a large family for which the method of [9] does not work ([10]).

In Section <ref> we present the classical `sum of two squares' problem in our 2D substitution context. This example also illustrates the fact that our technique extends from fixed points  of substitutions to 0-1-valued morphic sequences, i.e.,  fixed points of substitutions on arbitrary alphabets which are mapped to $\{0,1\}^*$ by a letter-to-letter map.

§ SUBSTITUTIONS AND THEIR PRODUCTS

A substitution $\sigma$ is a  homomorphism of the monoid $\A^*$ of words over an alphabet $\A$, that is, $\A^*$ consists of concatenations $a_1\dots a_m$ with $a_i$ from $\A$ for $i=1\dots m$, and $\sigma(vw)=\sigma(v)\sigma(w)$ for all $v,w\in \A^*$.
Since we are interested in the characteristic functions of 0-1-words, we simplify the presentation and take $\A=\{0,1\}$.

Let $\sigma$ on $\A^*=\{0,1\}^*$ be a substitution given by

$$\sigma(0)=a_1\dots a_m, \quad \sigma(1)=b_1\dots b_{m'},$$

for two natural numbers $m$ and $m'$. We then define the direct product substitution $\sigma\!\times\!\sigma$ on the alphabet $\A\times \A$ by

\begin{align*}
\sigma\!\times\!\sigma((0,0))_{(k,\ell)}  & =(a_k,a_\ell), \quad \sigma\!\times\!\sigma((0,1))_{(k,\ell')}=(a_k,b_{\ell'}),  \\
\sigma\!\times\!\sigma((1,0))_{(k',\ell)} & =(b_{k'},a_\ell),   \quad \sigma\!\times\!\sigma((1,1))_{(k',\ell')}=(b_{k'},b_{\ell'}),
\end{align*}

where $1\le k \le m, 1\le \ell \le m, 1\le k' \le m', 1\le \ell' \le m'.$

Example  Let $\sigma$ be the Fibonacci substitution given by $\sigma(0)=01,\; \sigma(1)=0$ then, abbreviating $(i,j)$ to $ij$ for $i,j=0,1$
the direct product substitution $\sigma\!\times\!\sigma$ is given  by

[scale=.3,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $00$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $01$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $10$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $11$
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P10);
#1=10(P00);  (P01);
#1=11(P00);  ;

(0,0) coordinate (a);   (a);  at (2,0) $\mapsto$;  (4,0) coordinate (t00);  θ(00,t00);
(8,0) coordinate (b);   (b);  at (10,0) $\mapsto$;  (12,0) coordinate (t01);  θ(01,t01);
(16,0) coordinate (c);  (c);  at (18,0) $\mapsto$;  (20,0) coordinate (t10); θ(10,t10);
(23,0) coordinate (d);  (d);  at (25,0) $\mapsto$; (27,0) coordinate (t11); θ(11,t11);

Suppose that $\sigma(0)$ has prefix 0. Then $\sigma^n(0)$ converges to a fixed point $x$ of $\sigma$ as $n\rightarrow\infty$, and so
$[\sigma\!\times\!\sigma]^n((0,0))$ converges to the 2D fixed point $x\times x$ of $\sigma\!\times\!\sigma$.

Thus the element of $x\times x$  with index $(k,\ell)$  contains the symbol $(0,0)$ if and only if $x_k=0, x_\ell=0$ if and only if  $k\in A$ and $\ell\in A$.

Let the diagonal words $d_n$ be defined for $n\ge 2$ by

\begin{equation}\label{eq:diag}
d_n =\{[x\times x]_{(k,\ell)}: k+\ell =n, \, k\ge 1, \, l\ge 1\}.
\end{equation}

Then $n\in A+A$ if and only if $d_n$ contains a symbol $(0,0)$. So we have obtained the following theorem.

Let $\sigma$ be a substitution on $\{0,1\}$, such that $\sigma(0)$ has prefix 0, and let $x$ be the fixed point of $\sigma$ with prefix $0$. Let $A$ be the sequence of positions of $0$ in $x$. Let $d_n$ be the diagonal words occurring in $x\times x$, for $n\ge 2$.
Then $A+A= \{n\ge 2:  \: the\: symbol\:  (0,0)\:  occurs\: in\: d_n\}.$

Note that we also have that if $B$ is the sequence of positions of $1$ in $x$, then  $B+B= \{n\ge 2:  \: the\: symbol\:  (1,1)\:  occurs\: in\: d_n\}$, and $A+B= \{n\ge 2:  \: the\: symbol\:  (0,1)\:or\: (1,0)\:  occurs\: in\: d_n\}$.

§ THE FIBONACCI WORD

The Fibonacci word $x_{\rm F}$ is the unique fixed point of the substitution $\sigma:\, 0\mapsto 01, \, 1\mapsto 0$.
Let $\varphi:=\sigma\!\times\!\sigma$ be the direct product of $\sigma$.

It is convenient to code   $a:=(0,0), b:=(0,1), c:=(1,0), d:=(1,1)$, and to assign colors to these four letters so that one obtains more easily an idea of the structure of $\varphi^n(a)$. The direct product substitution $\varphi$ on this new alphabet is given by

[scale=.3,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P10);
#1=10(P00);  (P01);
#1=11(P00);  ;

(0,0) coordinate (a);   (a);  at (2,0) $\mapsto$;  (4,0) coordinate (t00);  θ(00,t00);
(8,0) coordinate (b);   (b);  at (10,0) $\mapsto$;  (12,0) coordinate (t01);  θ(01,t01);
(16,0) coordinate (c);  (c);  at (18,0) $\mapsto$;  (20,0) coordinate (t10); θ(10,t10);
(23,0) coordinate (d);  (d);  at (25,0) $\mapsto$; (27,0) coordinate (t11); θ(11,t11);

Here are some examples of 2D words obtained by  iterating  $\varphi$.

[scale=.3,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P10);
#1=10(P00);  (P01);
#1=11(P00);  ;
(#1,#2)(#2) coordinate (Q0);
(Q0)++(0,0) coordinate  (Q00); (Q0)++(2,0) coordinate (Q10); (Q0)++(0,2) coordinate (Q01); (Q0)++(2,2) coordinate (Q11);
#1=00θ(00,Q00);  θ(01,Q01); θ(10,Q10); θ(11,Q11);
#1=01θ(00,Q00);  θ(10,Q10);
#1=10θ(00,Q00);  θ(01,Q01);
#1=11θ(00,Q00);  ;
(#1,#2)(#2) coordinate (R0);
(R0)++(0,0) coordinate  (R00); (R0)++(3,0) coordinate (R10); (R0)++(0,3) coordinate (R01); (R0)++(3,3) coordinate (R11);
#1=00(00,R00);  (01,R01); (10,R10); (11,R11);
#1=01(00,R00);  (10,R10);
#1=10(00,R00);  (01,R01);
#1=11(00,R00);  ;
(#1,#2)(#2) coordinate (S0);
(S0)++(0,0) coordinate  (S00); (S0)++(5,0) coordinate (S10); (S0)++(0,5) coordinate (S01); (S0)++(5,5) coordinate (S11);
#1=00(00,S00);  (01,S01); (10,S10); (11,S11);
#1=01(00,S00);  (10,S10);
#1=10(00,S00);  (01,S10);
#1=11(00,S00);  ;
(00,(0,0)); (00,(9,0)); (00,(20,0)); (01,(34,0)); (01,(43,0));
[circle, inner sep=1.4pt] at (-2.5,0) $\varphi^2(a)\!:$ ;
[circle, inner sep=1.4pt] at (6.5,0) $\varphi^3(a)\!:$ ;
[circle, inner sep=1.4pt] at (17.5,0) $\varphi^4(a)\!:$ ;
[circle, inner sep=1.4pt] at (31.5,0) $\varphi^2(b)\!:$ ;
[circle, inner sep=1.4pt] at (40.5,0) $\varphi^3(b)\!:$ ;

Note that

\begin{equation}\label{eq:diag}
d_2=a, \:d_3=bc, \:d_4=ada, \:d_5=acba,\: d_6=bcabc.
\end{equation}

We are now in a position to give a completely different proof of Theorem 3.1. of [7].

Let $A=1, 3, 4, 6, 8, 9, 11, 12, 14, 16,\dots$ be the  sequence $([n\alpha])$, where $\alpha$ is the golden mean. Then

Proof: Recall from the introduction that the sequence $([n\alpha])$ is equal to the sequence of positions of 0 in the infinite Fibonacci word $x_{\bf F}$.

We use Theorem <ref>. Trivially, 1 is not in $A+A$, and we see that $a$ does not occur in $d_3$. Further we see that $a$ does occur in $d_2$, $d_4$ and $d_5$. So it remains to show that all $d_n$ contain a symbol  $a$ for $n\ge 6$.
To see this, note that the first fife columns on the left border of $x_{\rm F}\times x_{\rm F}$ are a concatenation of the 2D words $\varphi^3(a)$ and $\varphi^3(b)$ for all $n$.
So these five columns are composed of blocks $\varphi^3(a)$ on top of $\varphi^3(a)$, $\varphi^3(a)$ on top of $\varphi^3(b)$, and $\varphi^3(b)$ on top of $\varphi^3(a)$:

[scale=.3,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P10);
#1=10(P00);  (P01);
#1=11(P00); ;
(#1,#2)(#2) coordinate (Q0);
(Q0)++(0,0) coordinate  (Q00); (Q0)++(2,0) coordinate (Q10); (Q0)++(0,2) coordinate (Q01); (Q0)++(2,2) coordinate (Q11);
#1=00θ(00,Q00);  θ(01,Q01); θ(10,Q10); θ(11,Q11);
#1=01θ(00,Q00);  θ(10,Q10);
#1=10θ(00,Q00);  θ(01,Q01);
#1=11θ(00,Q00);  ;
(#1,#2)(#2) coordinate (R0);
(R0)++(0,0) coordinate  (R00); (R0)++(3,0) coordinate (R10); (R0)++(0,3) coordinate (R01); (R0)++(3,3) coordinate (R11);
#1=00(00,R00);  (01,R01); (10,R10); (11,R11);
#1=01(00,R00);  (10,R10);
#1=10(00,R00);  (01,R01);
#1=11(00,R00); ;

[circle, inner sep=2.4pt] at (-23,62)  $\varphi^3(a):$; [circle, inner sep=2.4pt] at (-23,58)  $\varphi^3(a):$;
[circle, inner sep=2.4pt] at (-13,62)  $\varphi^3(a):$; [circle, inner sep=2.4pt] at (-13,58)  $\varphi^3(b):$;
[circle, inner sep=2.4pt] at (-3,63)  $\varphi^3(b):$; [circle, inner sep=2.4pt] at (-3,59)  $\varphi^3(a):$;
(-20,60)++(-1/2,-1/2) coordinate (left); (left)++(5,0) coordinate (right); [draw=gray,very thick] (left)–(right);
(-10,60)++(-1/2,-1/2) coordinate (left); (left)++(5,0) coordinate (right); [draw=gray,very thick] (left)–(right);
(-0,62)++(-1/2,-1/2) coordinate (left); (left)++(5,0) coordinate (right); [draw=gray,very thick] (left)–(right);

In all three cases, the diagonal words starting at the left border in the top $\varphi^3(a)$ or $\varphi^3(b)$ will cross a square with a symbol $a$ in one of the four left most columns, which finishes the proof. $\Box$.

Remark  The proof shows that if $n=a+a'$, with $a,a'$ from $A$, then  $a$ can always be chosen from the set $\{1,3,4\}$.

§ THE THUE-MORSE WORD

The Thue-Morse word $t$ is the fixed point of the substitution $\theta:\,0\mapsto 01, 1\mapsto 10$. Although it is general practice (as in [7]) to index Beatty sequences starting from $n=1$, and similarly for fixed points of substitutions, this is not the case for the Thue-Morse sequence.
The Thue-Morse word $t=t_0t_1\dots=01101001\dots$ is indexed starting from $n=0$. Thus  $A=0,3,5,6,\dots$ gives the positions of 0 in $t$, and $B=1,2,4,7,\dots$ the positions of 1 in $t$.

Let $\overline{0}=1, \overline{1}=0$ be the symmetry operator on $\{0,1\}^*$. Note that $\theta$ is symmetric, i.e., $\overline{\theta(i)}=\theta(\overline{i})$ for $i=0,1$.

The direct product substitution of $\theta$ is the 2D substitution $\tau$  defined on the symbols $(i,j)$ for $i,j=0,1$ (written as $ij$) by

\begin{align*}\label{eq:tau}
& i\overline{j}  \;\,  \overline{i}\overline{j}   \; &     \\[-.1cm]
\tau:\;ij\mapsto\; &  ij            \;\,  \overline{i}j.
\end{align*}

When we code    $  a:=00, b:=01, c:=10, d:=11$, and  color the squares of the symbols $a,b,c$ and $d$, then the substitution $\tau$ takes the form

[scale=.4,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P01); (P10); (P11);
#1=10(P00);  (P01); (P10); (P11);
#1=11(P00);  (P01); (P10); (P11);
(0,0) coordinate (ZZ);   (ZZ);  at (1.5,0) $\mapsto$;  (3,0) coordinate (t00);  θ(00,t00);
(7,0) coordinate (ZO);   (ZO);  at (8.5,0) $\mapsto$;  (10,0) coordinate (t01);  θ(01,t01);
(14,0) coordinate (OZ);  (OZ);  at (15.5,0) $\mapsto$;  (17,0) coordinate (t10); θ(10,t10);
(21,0) coordinate (OO);  (OO);  at (22.5,0) $\mapsto$; (24,0) coordinate (t11); θ(11,t11);
[circle, inner sep=1.4pt] at (5,0)  ,; [circle, inner sep=1.4pt] at (12,0)  ,;
[circle, inner sep=1.4pt] at (19,0)  ,; [circle, inner sep=1.4pt] at (26,0)  .;

Since we start $t$ at index 0, we have to reconsider the definition of the diagonal words. It turns out that it is convenient to index these by their lengths. So $d_1=a, d_2=bc$, etc.
The iterate $\tau^8(a)$ with the diagonal words $d_1,d_2,\dots,d_{16}$ indicated by lines, is given by

[scale=.4,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$

θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P01); (P10); (P11);
#1=10(P00);  (P01); (P10); (P11);
#1=11(P00);  (P01); (P10); (P11);

(#1,#2)(#2) coordinate (Q0);
(Q0)++(0,0) coordinate  (Q00); (Q0)++(2,0) coordinate (Q10); (Q0)++(0,2) coordinate (Q01); (Q0)++(2,2) coordinate (Q11);
#1=00θ(00,Q00);  θ(01,Q01); θ(10,Q10); θ(11,Q11);
#1=01θ(01,Q00);  θ(00,Q01); θ(11,Q10); θ(10,Q11);
#1=10θ(10,Q00);  θ(11,Q01); θ(00,Q10); θ(01,Q11);
#1=11θ(11,Q00);  θ(10,Q01); θ(01,Q10); θ(00,Q11);

(#1,#2)(#2) coordinate (R0);
(R0)++(0,0) coordinate  (R00); (R0)++(4,0) coordinate (R10); (R0)++(0,4) coordinate (R01); (R0)++(4,4) coordinate (R11);
#1=00(00,R00);  (01,R01); (10,R10); (11,R11);
#1=01(01,R00);  (00,R01); (11,R10); (10,R11);
#1=10(10,R00);  (11,R01); (00,R10); (01,R11);
#1=11(11,R00);  (10,R01); (01,R10); (00,R11);

(#1,#2)(#2) coordinate (S0);
(S0)++(0,0) coordinate  (S00); (S0)++(8,0) coordinate (S10); (S0)++(0,8) coordinate (S01); (S0)++(8,8) coordinate (S11);
#1=00(00,S00);  (01,S01); (10,S10); (11,S11);
#1=01(01,S00);  (00,S01); (11,S10); (10,S11);
#1=10(10,S00);  (11,S01); (00,S10); (01,S11);
#1=11(11,S00);  (10,S01); (01,S10); (00,S11);
(#1) coordinate (start);
ıin 1,2,...,11 (start)++(12,0) coordinate (start);;;
(#1) coordinate (start);
ıin 1,2,...,16 (start)++(0,ı) coordinate (up); (start)++(ı,0) coordinate (down); [gray] (up)–(down); ;

Our goal is to describe the diagonal words $d_n$, for $n=1,2\dots$. We see that in particular

\begin{equation}\label{eq:diagd15}
d_{1}=a, \:d_2=bc, \:d_3=bdc, \:d_4=adda.
\end{equation}

Let $\sigma$ be the morphism given by $\sigma(ij) = i\overline{j}, \overline{i}j$ for $i,j=0,1$. Then

$$d_{2n}=\sigma(d_n), \quad  for\; n=1,2,\dots.$$

Let $\beta$ be the 2-to-1 morphism given by $\beta(ij,i'j')=ij'$ for $i,j,i',j'=0,1$. Then

$$d_{2n+1}=\beta(d_{2n+2}), \quad  for\; n=1,2,\dots.$$

Proof: Let the lines $\Delta_n:=\{(x,y):x+y=n-1\}$ correspond to the diagonal words $d_n$. The line $\Delta_{2n}$ cuts through the diagonals of the $2\times 2$-blocks $\tau(a),\dots,\tau(d)$, and these $2\times 2$-blocks are images under $\tau$ of the symbols $a,\dots,d$ on the line $\Delta_n$. Thus  $d_{2n}=\sigma(d_n)$ now follows from the observation that the $i\overline{j}, \overline{i}j$ are exactly the symbols on the diagonals of the $2\times 2$ blocks $\tau(ij)$.

For the second statement, we look at the successive $2\times 2$ blocks $\tau(a),\dots,\tau(d)$  cut by  the line $\Delta_{2n+2}$. The line $\Delta_{2n+1}$ cuts through the lower left squares of these blocks. This implies that the symbols at the odd indices of the word $d_{2n+1}$  can be read off directly from the successive pairs of symbols of  $d_{2n+2}$ by mapping the diagonals of the $2\times 2$ $\tau$-blocks to the symbols at the lower left corner of these blocks.
This gives   $\beta(i\overline{j},\,\overline{i}j)=   ij$, which is an instance of $\beta(ij,i'j')=ij'$.

For the symbols at the even indices of the word $d_{2n+1}$ we have to do more work: we have to determine the symbol at the upper right corner of the $\tau$-block directly below a pair of $\tau$-blocks from the pair of symbols situated at lower left and the upper right corner of these blocks. The situation is as in the following drawing:

[scale=.5,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); [draw=black]  (shP0) rectangle +(1,1);
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); [draw=black]  (shP0) rectangle +(1,1);
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); [draw=black]  (shP0) rectangle +(1,1);
(#1)(#1)+(-1/2,-1/2) coordinate (shP0); [draw=black]  (shP0) rectangle +(1,1);
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P01); (P10); (P11);
#1=10(P00);  (P01); (P10); (P11);
#1=11(P00);  (P01); (P10); (P11);
(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00); (P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11);
#1=10(P00);  (P01); (P10); (P11);

(6,4) coordinate (OZ);  (10,OZ);  (8,2) coordinate (ZZ); (10,ZZ);  (6,2) coordinate (OZ); (10,OZ);
[circle, inner sep=1.4pt] at (8,3)  $i'\overline{j'}$;
[circle, inner sep=5.4pt] at (7,4)  $\overline{i}j$;

(12,4) coordinate (OZ);  (10,OZ);  (14,2) coordinate (ZZ); (10,ZZ);  (12,2) coordinate (OZ); (10,OZ);
[circle, inner sep=1.4pt] at (14,3)  $i'\overline{j'}$ ;
[circle, inner sep=5.4pt] at (13,4)  $\overline{i}j$;
[circle, inner sep=5.4pt] at (13,3)  $\overline{i}\,\overline{j'}$;
[circle, inner sep=5.4pt] at (10.5,3)  $\Rightarrow$;

This gives: $\beta(\overline{i}j, \footnotesize i'\overline{j'}) = \overline{i}\,\overline{j'}$,
which happens to be a version of $\beta(ij,i'j')=ij'$.     $\Box$

Let $A=0,3,5,6,9,10,\dots$ give the positions of $0$'s in Thue Morse sequence, and $B=1,2,4,6,7,8,\dots$ the positions of $1$ in the Thue Morse sequence. Then

$$A+A=\mathbb{N}_0\setminus\{2,4,\,2^{2n+1}\!-1, n\ge 0\}.$$

$$B+B=\mathbb{N}_0\setminus\{2^{2n+1}\!-1, n\ge 0\}.$$

$$A+B=\mathbb{N}_0\setminus\{2^{2n}\!-1, n\ge 0\}.$$

Proof: We only prove the $A+A$ result. The other two can be proved in a similar way.

According to Theorem <ref>,   $n\in A+A$ if and only if $d_{n-1}$ contains a symbol $(0,0)=a$. Note that we had to shift the diagonal words, as they are redefined for this section.

Let $\sigma$ be the morphism from Proposition <ref> in $abcd$-coding,  $\sigma(a)=bc, \,\sigma(b)=ad, \,\sigma(c)=da, \,\sigma(d)=cb$. Then $\sigma^2$ is given by

$$\sigma^2(a)=adda, \,\sigma^2(b)=bccb, \,\sigma^2(c)=cbbc, \,\sigma^2(d)=daad.$$

It follows from this and Equation (<ref>) that $d(2^{2n+1})\in \{b,c\}^*,  \,d(2^{2n})\in \{a,d\}^*$, for $n\ge 0$.

This takes care of the $2^{2n+1}\!-1$ part of the theorem.  Obviously 2 and 4 are not in $A+A$. It remains to prove that all other numbers are in $A+A$, which is equivalent to proving that a symbol $a$  occurs in all  diagonal words  $d_n$ with $n>5$ and $n$ not equal to a power of 2.
It looks attractive to use Proposition <ref> to accomplish this. Indeed, let $T$ be the set consisting of the twelve 2-blocks $ab, ac,\dots,dc$ with two different symbols. Then one may check that all elements of $T$ occur in the length 4 blocks from $\sigma(T)$, and also in the length 3 blocks from $\beta(\sigma(T))$. This implies that a diagonal word $d_n$ in which all elements from $T$ occur, propagates this property to $d_{2n}$ and $d_{2n-1}$. However, there are many diagonal words in which not all blocks from $T$ occur, in fact all $n$ of the form $n=2^N+2^K$ for some non-negative integers $K$ and $N$. This makes an induction proof very complex, and so we will give a completely different proof.

In Figure <ref> we depict the  $\tau^3$-squares. Here the red lines indicate diagonal words without   the symbol $a$, and the  green lines indicate diagonal words  with a symbol $a$.

Diagonals in $\tau^3$-squares.

We also indicated (parts of) the diagonal words above the main diagonal. These are denoted by  $d_1^+,\dots,d_7^+$, starting from the main diagonal. For instance the block $\tau^3(a)$ has the two red diagonal lines $d_6^+$ and $d_7^+$.

Diagonals in $\tau^4$-squares inherited from $\tau^3$-blocks.

Figure <ref> gives the red and green line structure for the $\tau^4$-blocks inherited from the $\tau^3$-blocks.

Final color assignment to the diagonals in $\tau^4$-blocks.

Observe that any diagonal that is partly red, partly green, actually should be a
green diagonal. This leads to Figure <ref>.

We finish the proof with induction. Starting from $n=4$, the blocks $\tau^n(a), \tau^n(b), \tau^n(c), \tau^n(d)$ have the following properties:

1) There are some red diagonals $d_i$ with $i\in \{1,2,3,4,5\}$,

2) There are some red diagonals $d^+_i$ with $i\in \{2^n-5,\dots,2^n-1\}$,

3) There are some red diagonals $d_i$ and $d^+_i$ with $i\in \{2^N, N\ge 2 \}$,

4) All other diagonals $d_i$ and $d^+_i$ are green.

One checks that these four properties hold for $n=4$, and then using the induction hypothesis makes the step from $n$ to $n+1$ in the same way as the step from $n=3$ to $n=4$ is made as in Figures <ref>, <ref> and <ref>. $\Box$.

§ THE VON NEUMANN WORD

The von Neumann word $u$ was introduced in the paper [1]. One has $u=u_0u_1\dots=1101100110110\dots$, fixed point of the substitution  $\nu:\quad 0\mapsto 0,\; 1\mapsto 110.$

The 2D von Neumann substitution $\psi=\nu\times\nu$ is  given by

[scale=.4,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $a$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $b$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $c$
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!30!white, draw=black]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $d$
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate (P00); (P0)++(0,1) coordinate (P01); (P0)++(0,2) coordinate (P02);
(P0)++(1,0) coordinate (P10); (P0)++(1,1) coordinate (P11); (P0)++(1,2) coordinate (P12);
(P0)++(2,0) coordinate (P20); (P0)++(2,1) coordinate (P21); (P0)++(2,2) coordinate (P22);
#1=00(P00) ;
#1=01(P00); (P01); (P02) ;
#1=10(P00); (P10); (P20) ;
#1=11(P00); (P01); (P02); (P10); (P11); (P12); (P20); (P21); (P22) ;
(0,-6) coordinate (a);   (a);  at (2,-6) $\mapsto$;  (4,-6) coordinate (t00);   θ(00,t00);
(7,-6) coordinate (b);   (b);  at (9,-6) $\mapsto$;  (11,-6) coordinate (t01);  θ(01,t01);
(14,-6) coordinate (c);  (c);  at (16,-6) $\mapsto$; (18,-6) coordinate (t10);  θ(10,t10);
(23,-6) coordinate (d);  (d);  at (25,-6) $\mapsto$; (27,-6) coordinate (t11);  θ(11,t11);

Before we continue, we formulate a simple lemma on the language $L_\nu$ of the substitution $\nu$, i.e., on the collection of words that can occur in $u$.

The words $010, 111$ and $101101$ do not occur in $L_\nu$.

Using Lemma <ref> one simply derives the next lemma.

The word $11$ occurs uniquely as suffix of $011$. The word $101$ occurs uniquely as suffix of $001101$. The word $10110$ occurs uniquely as suffix of $0110110$.

The prefix tree for the von Neumann words.

Diagonals for 2D von Neumann words in Case 1,2,3.

We split the possible occurrences of words of $L_\nu$ in the von Neumann word $u$ according to a tree of prefixes, dividing these in six possible cases. This prefix tree is given in Figure <ref>, where $\Lambda$ is the empty word. Note that we used Lemma <ref> again to determine the offspring of the node labeled 10.

Let $A=3,6,7,10,13,\dots$ give the positions of $0$'s in the von Neumann sequence, and $B=1,2,4,5,8,9\dots$ the positions of $1$ in the von Neumann sequence. Then

Proof: The proof of $B+B=\mathbb{N}\setminus\{1\}$ is almost trivial.

The proof for $A+A$ is based on Theorem <ref>, with the trivial change that we consider the fixed point starting with 1. The proof runs as the proof for the Fibonacci word, except that here we need the 15 left most columns.

First one checks easily that the numbers $2,3,4,5,7,8,11,15$ are not in $A+A$.

To handle the diagonal words $d_n$ for $n\ge 16$, we split into the six cases we introduced in Figure <ref>. In Case 2, 3 and 4 we consider the three diagonals starting at the left border in the blocks $\psi(d)$, in the remaining four the single diagonal starting in $\psi(c)$. See Figure <ref> for the three cases 00, 11 and 001. This means that in Case 1 the  top two blocks are $\psi(c)\psi(c)$, in Case 2 they are $\psi(d)\psi(d)$, and in Case 3 they are $\psi(c)\psi(c)\psi(d)$.

Diagonals for 2D von Neumann words in  Case 4,5,6.

In Case 2  a block $\psi(c)$, must precede the top blocks $\psi(d)\psi(d)$, by Lemma 2. In a similar way there are forced blocks in Case 4 and Case 6.

The other blocks are then inserted following the von Neumannn words $ccaccaa\dots$ or $ddbddbb\dots$. In all cases, except the last, a row $ccaccaacca\dots$ has been added at the bottom. This is allowed because both $\psi(c)$ and $\psi(d)$ have the letter $c$ as prefix.

The red diagonals are those parts of the three diagonals starting in the top block $\psi(d)$, respectively the single diagonal starting in $\psi(c)$,  that do not cross a square with label $a$, i.e., the red diagonal ends just before the first $a$-square.
See Figure <ref> for the three cases remaining cases.

The fact that in all six cases a square with label $a$ is encountered finishes the proof. $\Box$

§ THE SUM OF SQUARES

Here we give a classical example of a sumset. Let $A=\{n^2: n\ge 1\}$. It is well known (see, e.g., [3]), that the characteristic function of $A$, as a word, is a letter to letter substitution $\lambda$ of the fixed point with prefix 0 of the morphism

$$ 0\mapsto 01, \quad  1\mapsto 221, \quad  2\mapsto 2.$$

The letter to letter map is given by $\lambda(0)=0, \;\lambda(1)=1,\;\lambda(2)=0$.

The corresponding 2D substitution is the morphism $\mu$ given by

[scale=.4,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!40!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $00$;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!60!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $01$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $02$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!60!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $10$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!60!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $11$;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $12$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $20$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $21$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!10!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $22$ ;
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00);
(P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11); (P0)++(2,0) coordinate (P20); (P0)++(0,2) coordinate (P02); (P0)++(1,2) coordinate (P12); (P0)++(2,1) coordinate (P21); (P0)++(2,2) coordinate (P22);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P01); (P10); (P11); (P02); (P12);
#1=02(P00);  (P10) ;
#1=10(P00);  (P01); (P11); (P10);  (P21); (P20);
#1=11(P00);   (P01);  (P02);
(P10);   (P11);  (P12);
(P20);   (P21);  (P22) ;
#1=12(P00);  (P10); (P20);
#1=20(P00);  (P01); ;
#1=21(P00);  (P01); (P02) ;
#1=22(P00); ;

(0,-4) coordinate (ZZ);   (ZZ);  at (2,-4) $\mapsto$;  (4,-4) coordinate (t00);  θ(00,t00);
(10,-4) coordinate (ZO);   (ZO);  at (12,-4) $\mapsto$;  (14,-4) coordinate (t01);  θ(01,t01);
(20,-4) coordinate (ZT);  (ZT);  at (22,-4) $\mapsto$;  (24,-4) coordinate (t02);  θ(02,t02);
(0,-8) coordinate (OZ);  (OZ);  at (2,-8) $\mapsto$;  (4,-8) coordinate (t10); θ(10,t10);
(10,-8) coordinate (OO);  (OO);  at (12,-8) $\mapsto$; (14,-8) coordinate (t11); θ(11,t11);
(20,-8) coordinate (OT);  (OT);  at (22,-8) $\mapsto$; (24,-8) coordinate (t12); θ(12,t12);
(0,-12) coordinate (TZ);  (TZ);  at (2,-12) $\mapsto$; (4,-12) coordinate (t20); θ(20,t20);
(10,-12) coordinate (TO);  (TO);  at (12,-12) $\mapsto$; (14,-12) coordinate (t21); θ(21,t21);
(20,-12) coordinate (TT);  (TT);  at (22,-12) $\mapsto$; (24,-12) coordinate (t22); θ(22,t22);

The numbers $s=n^2+m^2$ that are a sum of two squares occur as sums of the indices $(n^2,m^2)$  of the squares with the symbols $11$ in the fixed point $S$ of $\mu$, limit of $\mu^n(00)$ as $n\rightarrow \infty$.

At first sight it is surprising that this set has a fractal structure, but due
to the fact that the morphism $\mu$ is not primitive (i.e., the incidence matrix of the substitution is reducible), there is no exponential scaling structure in the 2D word $S$, but rather a polynomial one. The latter simply amounts to the recursion $(n+1)^2=n^2+(2n+1)$.
This recursion is clearly visible in  Figure <ref>, which displays the two-dimensional word  $\mu^4(00)$.

[scale=.4,rounded corners=0]
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=green!40!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $00$;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!60!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $01$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $02$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!60!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $10$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=red!60!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $11$;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $12$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=yellow!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $20$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!30!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $21$ ;
(#1)(#1)+(-1/2,-1/2) coordinate (shP0);
[fill=blue!10!white, draw=gray,thick]  (shP0) rectangle +(1,1); [circle, inner sep=5.4pt] at (#1)  $22$ ;
θ(#1,#2)(#2) coordinate (P0);
(P0)++(0,0) coordinate  (P00);
(P0)++(1,0) coordinate (P10); (P0)++(0,1) coordinate (P01); (P0)++(1,1) coordinate (P11); (P0)++(2,0) coordinate (P20); (P0)++(0,2) coordinate (P02); (P0)++(1,2) coordinate (P12); (P0)++(2,1) coordinate (P21); (P0)++(2,2) coordinate (P22);
#1=00(P00);  (P01); (P10); (P11);
#1=01(P00);  (P01); (P10); (P11); (P02); (P12);
#1=02(P00);  (P10) ;
#1=10(P00);  (P01); (P11); (P10);  (P21); (P20);
#1=11(P00);   (P01);  (P02);
(P10);   (P11);  (P12);
(P20);   (P21);  (P22) ;
#1=12(P00);  (P10); (P20);
#1=20(P00);  (P01); ;
#1=21(P00);  (P01); (P02) ;
#1=22(P00); ;
(#1,#2)(#2) coordinate (Q0);
(Q0)++(0,0) coordinate (Q00); (Q0)++(2,0) coordinate (Q10); (Q0)++(4,0) coordinate (Q20);
(Q0)++(0,1) coordinate (Q01); (Q0)++(2,1) coordinate (Q11); (Q0)++(4,4) coordinate (Q21);
(Q0)++(0,2) coordinate (Q02); (Q0)++(2,2) coordinate (Q12); (Q0)++(2,2) coordinate (Q22);
#1=00 (Q0)++(0,2) coordinate (Q01); (Q0)++(2,2) coordinate (Q11);
θ(00,Q00);  θ(01,Q01); θ(10,Q10); θ(11,Q11);
#1=01θ(02,Q00);  θ(02,Q01); θ(12,Q10); θ(12,Q11); θ(01,Q02); θ(11,Q12);
#1=02θ(02,Q00);  θ(12,Q10) ;
#1=10(Q0)++(1,0) coordinate (Q10); (Q0)++(2,0) coordinate (Q20);
(Q0)++(0,2) coordinate (Q01); (Q0)++(1,2) coordinate (Q11); (Q0)++(2,2) coordinate (Q21);
θ(20,Q00);  θ(20,Q10); θ(21,Q11); θ(21,Q01);  θ(11,Q21); θ(10,Q20);
#1=11(Q0)++(1,0) coordinate (Q10); (Q0)++(2,0) coordinate (Q20);
(Q0)++(0,1) coordinate (Q01); (Q0)++(1,1) coordinate (Q11); (Q0)++(2,1) coordinate (Q21);
(Q0)++(0,2) coordinate (Q02); (Q0)++(1,2) coordinate (Q12); (Q0)++(2,2) coordinate (Q22);
θ(22,Q00);  θ(22,Q01);  θ(21,Q02);
θ(22,Q10);  θ(22,Q11);  θ(21,Q12);
θ(12,Q20);  θ(12,Q21);  θ(11,Q22) ;
#1=12(Q0)++(1,0) coordinate (Q10); (Q0)++(2,0) coordinate (Q20);
θ(22,Q00);  θ(22,Q10); θ(12,Q20);
#1=20 (Q0)++(0,2) coordinate (Q01); θ(20,Q00);  θ(21,Q01); ;
#1=21θ(22,Q00);  θ(22,Q01); θ(21,Q02) ;
(#1,#2)(#2) coordinate (R0);
(R0)++(0,0) coordinate  (R00); (R0)++(5,0) coordinate (R10); (R0)++(0,5) coordinate (R01); (R0)++(5,5) coordinate (R11);
(R0)++(0,5) coordinate (R02); (R0)++(5,5) coordinate (R12); (R0)++(5,5) coordinate (R22);
#1=00(00,R00);  (01,R01); (10,R10); (11,R11);
#1=01(R0)++(5,0) coordinate (R10); (R0)++(0,1) coordinate (R01);
(R0)++(5,1) coordinate (R11); (R0)++(0,2) coordinate (R02);  (R0)++(5,2) coordinate (R12);
(02,R00);  (02,R01); (12,R10); (12,R11); (01,R02); (11,R12);
#1=02(02,R00);  (12,R10) ;
#1=10(R0)++(1,0) coordinate (R10); (R0)++(2,0) coordinate (R20);
(R0)++(0,5) coordinate (R01); (R0)++(1,5) coordinate (R11); (R0)++(2,5) coordinate (R21);
(20,R00);  (20,R10); (21,R11); (21,R01);  (11,R21); (10,R20);
#1=11(R0)++(1,0) coordinate (R10); (R0)++(2,0) coordinate (R20);
(R0)++(0,1) coordinate (R01); (R0)++(1,1) coordinate (R11); (R0)++(2,1) coordinate (R21);
(R0)++(0,2) coordinate (R02); (R0)++(1,2) coordinate (R12); (R0)++(2,2) coordinate (R22);
(22,R00);  (22,R01);  (21,R02);
(22,R10);  (22,R11);  (21,R12);
(12,R20);  (12,R21);  (11,R22) ;
#1=12(R0)++(1,0) coordinate (R10); (R0)++(2,0) coordinate (R20);
(22,R00);  (22,R10); (12,R20);
#1=20(20,R00); (21,R01); ;
#1=21(R0)++(0,1) coordinate (R01); (R0)++(0,2) coordinate (R02);
(22,R00);  (22,R01); (21,R02) ;
(#1,#2)(#2) coordinate (S0);
(S0)++(0,0) coordinate  (S00); (S0)++(10,0) coordinate (S10); (S0)++(0,10) coordinate (S01); (S0)++(10,10) coordinate (S11);
#1=00(00,S00);  (01,S01); (10,S10); (11,S11) ;

The 2D word $\mu^4(00)$.

§ ACKNOWLEDGEMENT

I am grateful to Jean-Paul Allouche and Jeffrey Shallit for their useful comments on this paper.

[1]  J.-P. Allouche, J. Betrema and J. O. Shallit,    Sur des points fixes de morphismes d'un monoïde libre,
RAIRO - Theoretical Informatics and Applications - Informatique Théorique et Applications  23 (1989), 235–249.

[3] A. Cobham, Uniform tag sequences, Math. Systems Theory 6 (1972), 164–192.

[4] F.M. Dekking, K. Simon and B. Székely. The algebraic difference  of two random Cantor sets: the Larsson
family,  Annals of Probability 39 (2011), 549–586.

[6] N. Priebe Frank, A primer on substitution tilings of the Euclidean plane, Expo. Math. 26 (2008), 295–326.

[7]  S. Kawsumarng, T. Khemaratchatakumthorn, P. Noppakaew and P. Pongsriiam,  Sumsets associated with Wythoff sequences and Fibonacci numbers, Periodica Mathematica Hungarica 82 (2021), 98–113.

[8] M. Lothaire,  Algebraic combinatorics on words, Encyclopedia of Mathematics and its Applications 90, Cambridge University Press, 2002.

[9] J. Shallit,  Sumsets of Wythoff Sequences, Fibonacci Representation, and Beyond, arXiv:2006.04177v1 [math.CO], 7 Jun 2020 (2020).

[10] J. Shallit.  Personal communication. March 2021.

[11] A. Rajasekaran, J. Shallit, and T. Smith. Additive number theory via
automata theory. Theor. Comput. Systems 64 (2020), 542–567.

[12] E. Zeckendorf, Réprésentation des nombres naturels par une somme de
nombres de Fibonacci ou de nombres de Lucas,  Bull. Soc. Roy. Sci. Liège 41 (1972), 179–182.

AMS Classification Numbers: 11B13, 68R10.


# Boundary Control for Suppressing Chaotic Response to Dynamic Hydrogen
Blending in a Gas Pipeline

Luke S. Baker  Anatoly Zlotnik  Applied Mathematics & Plasma Physics Group,
Los Alamos National Laboratory, Los Alamos, NM, 87545 (e-mail: {lsbaker,
azlotnik}@ lanl.gov).

###### Abstract

It is known that periodic forcing of nonlinear flows can result in a chaotic
response under certain conditions. Such non-periodic and chaotic solutions
have been observed in simulations of heterogeneous gas flow in a pipeline with
periodic, time-varying boundary conditions. In this paper, we examine a
proportional feedback law for boundary control of a parabolic partial
differential equation system that represents the flow of two gases through a
pipe. We demonstrate that periodic variation of the mass fraction of the
lighter gas at the pipe inlet can result in the chaotic propagation of gas
pressure waves, and show that appropriate flow control can suppress this
response. We examine phase space solutions for the single pipe system subject
to boundary control, and use numerical experiments to characterize conditions
for the controller gain to suppress chaos.

###### keywords:

Control of bifurcation and chaos, Control of renewable energy resources,
Energy control in transportation, Output feedback control, Control of fluid
flows and fluids-structures interactions

††thanks: This study was funded by the LANL Laboratory Directed R&D Project
“Efficient Multi-scale Modeling of Clean Hydrogen Blending in Large Natural
Gas Pipelines to Reduce Carbon Emissions”. Research conducted at Los Alamos
National Laboratory is done under the auspices of the National Nuclear
Security Administration of the U.S. Department of Energy under Contract No.
89233218CNA000001.

##  1 Introduction

Chaotic behavior has been observed in a variety of continuum systems that
involve gas dynamics, which are defined by partial differential equation (PDE)
systems. Prominent examples in nature include seasonal resonances and weather
patterns in the oceanic atmosphere related to cyclical wind bursts (Tziperman
et al. (1994); Eisenman et al. (2005)). Examples in engineered systems include
flame combustion of mixtures of hydrogen and air in microchannels (Pizza et
al. (2008); Alipoor and Mazaheri (2016)). The chaotic dynamics in these
systems are observed as high sensitivity to small perturbations in initial
conditions that result in rapid divergence of trajectories, and as irregular
responses to periodic boundary conditions.

Although chaos can be directly exploited in certain applications, such as
cryptography (Kocarev (2001)), it is in general not desirable in most
engineering settings, and in particular those that involve fluid flows. The
suppression of chaotic dynamics is desired in engineered systems with
observable states and controllable inputs, in order to maintain predictable
function within design limits. The application of control methods to suppress
chaos in a dynamical system has been examined in various studies (Boccaletti
et al. (2000); Zhang et al. (2009)). Approaches including open-loop control
(Dudnik et al. (1983)), feedback control (Pyragas (1992); Hu et al. (1995);
Schuster and Stemmler (1997)), and the Ott-Grebogi-Yorke (OGY) method (Ott et
al. (1990)) have all demonstrated the ability to control or suppress chaotic
behavior in dynamic systems.

Recently, the emergence of chaotic flows was shown to occur in the downstream
pressure of a pipeline that transports highly heterogeneous mixtures of gases
(Baker et al. (2023b)). The study was motivated from previous observations
that a sudden injection of hydrogen into a natural gas pipeline could cause
pressure, density, and mass flux values to diverge significantly from their
nominal natural gas values (Melaina et al. (2013); Eames et al. (2022);
Zlotnik et al. (2023)). Modeling flow phenomena for heterogeneous gas mixtures
is considerably more complex than for a single gas. The flow of natural gas in
a pipe is governed by a state equation and two PDEs that describe mass and
momentum conservation. Another PDE must be added to track the concentration of
each additional gas constituent. Because societies worldwide are investing in
clean energy resources for a transition to zero carbon economies, and hydrogen
is a compelling energy carrier that can be produced using renewable
electricity (Götz et al. (2016); Ozturk and Dincer (2021)), hydrogen blending
into gas pipelines is being evaluated to decarbonize energy systems while
using capital investments in infrastructure for their planned lifetimes (Hafsi
et al. (2019); Subani et al. (2017); Chaczykowski et al. (2018); Elaoud et al.
(2017)). While model predictive optimal control methods for pipeline
operations have been extended to the setting of hydrogen blending (Baker et
al. (2023a)), the possibility of chaotic pressure waves would challenge
contemporary pipeline operation methods (Baker et al. (2023b)), which use gas
compressors and pressure regulators to maintain predictable flows. This
compels methods to characterize the onset of chaotic dynamics in pipeline
systems, and control techniques to suppress such responses.

In this study, we review techniques for characterizing regions in the periodic
forcing parameter space that separate chaotic and non-chaotic responses to
variations of hydrogen mass fraction in the gas mixture entering a pipeline.
Sinusoidal variation of hydrogen concentration at the boundary of a pipeline
can create chaotic flows with large spike surges in pressure, density, mass
flux, and concentration variables that propagate along the pipe, and this was
demonstrated according to the rigorous mathematical definition of chaos Baker
et al. (2023b). Here, we apply proportional feedback control to compressor
settings in order to regulate the mass flow at the pipe boundary, and
demonstrate with numerical simulations that this controller can suppress
chaotic dynamics along the pipe.

The rest of this paper is organized as follows. In Section 2, the PDE system
that governs the flow of a mixture of gases in a pipe is introduced and
reduced to a finite-dimensional ordinary differential equation (ODE) system by
applying Chebyshev discretization in space. The method that we use to quantify
chaotic dynamics in the reduced ODE system is described in Section 3. In
Section 4, we define the proportional feedback control formula, in addition to
mathematically defining the suppression of chaotic behavior. Finally, Section
4 describes numerical experiments in which we demonstrate the performance of
the controller to suppress chaos over the space of periodic forcing
parameters.

##  2 Heterogeneous Gas Flow in a Pipe

The time-dependent flow of a mixture of natural gas and hydrogen in a
horizontal pipe is modeled with the parabolic PDE system (Chaczykowski et al.
(2018); Baker et al. (2023b))

$\displaystyle\frac{\partial\rho^{(m)}}{\partial t}+\frac{\partial}{\partial
x}\left(\frac{\rho^{(m)}}{\rho^{(1)}+\rho^{(2)}}\varphi\right)$
$\displaystyle=$ $\displaystyle 0,$ (1) $\displaystyle\frac{\partial}{\partial
x}\left(\sigma^{2}_{1}\rho^{(1)}+\sigma^{2}_{2}\rho^{(2)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda}{2D}\frac{\varphi|\varphi|}{\rho^{(1)}+\rho^{(2)}},\quad$
(2)

where (1) is defined for each constituent gas ($m=1$ for natural gas and $m=2$
for hydrogen). The variables are natural gas density $\rho^{(1)}(t,x)$,
hydrogen density $\rho^{(2)}(t,x)$, and total mass flux $\varphi(t,x)$ for
$t\in[0,T]$ and $x\in[0,\ell]$, where $T$ denotes the time horizon and $\ell$
denotes the length of the pipe. The diameter and friction factor of the pipe
are denoted by $D$ and $\lambda$. The wave speeds in natural gas and hydrogen
are denoted by $\sigma_{1}$ and $\sigma_{2}$, respectively. From Dalton’s law,
the total pressure $p$ of the mixture of hydrogen and natural gas is equal to
the summation of their partial pressures. For an ideal equation of state, the
total pressure is given by
$p=\sigma_{1}^{2}\rho^{(1)}+\sigma_{2}^{2}\rho^{(2)}$. The dynamics in (2)
assume this ideal relation. Because the ideal equation of state qualitatively
reflects the dynamics seen in the non-ideal case, we would expect to observe
phenomenologically similar flow behavior with the use of more accurate
nonlinear equations of state (Botros and Jensen (2022)).

Boundary conditions for equations (1)-(2) are defined by

$\displaystyle\rho^{(m)}(t,0)$ $\displaystyle=$ $\displaystyle
u(t)s^{(m)}(t),\qquad(m=1,\;2),$ (3) $\displaystyle\varphi(t,\ell)$
$\displaystyle=$ $\displaystyle\overline{q},$ (4)

where $s^{(m)}$ ($m=1$ and $m=2$) are known functions that represent the
partial densities of the respective gases at the inlet source, and
$\overline{q}$ is a known constant value of the total mass flux of the mixture
being withdrawn from the outlet of the pipe. The function $u$ is the boundary
control variable that represents the pressure increase after flow through a
compressor station located at the pipe inlet. To more clearly distinguish
between chaos and usual transient flow behavior, we specify initial conditions
and all boundary parameters, except for concentration, to be as steady as
possible. Thus we choose the source pressure
$\overline{p}=\sigma_{1}^{2}s^{(1)}+\sigma_{1}^{2}s^{(1)}$ and the withdrawal
flux $\overline{q}$ to be constant values in time. Moreover, the initial
condition is chosen to be the steady-state solution. As such, time-
dependencies of the flow variables result only from forcing the concentration
to fluctuate at the boundary. Hydrogen concentration (mass fraction) at the
boundary is defined by $\gamma=s^{(2)}/(s^{(1)}+s^{(2)})$. If $\gamma$ and
$\overline{p}$ are specified, then we may use their defining equations to
determine necessary $s^{(m)}$ that would achieve such an outcome. Henceforth,
we assume that $\overline{p}$ and $\gamma$ are specified, and we define

$\displaystyle s^{(1)}(t)=(1-\gamma(t))\overline{p}/\sigma_{1}^{2},\qquad
s^{(2)}(t)=\gamma(t)\overline{p}/\sigma_{2}^{2}.$ (5)

These definitions are enforced in the boundary conditions.

In general, the solution of (1)-(4) must be obtained numerically. We apply a
Chebyshev discretization scheme over the space domain $[0,\ell]$. The interval
$[0,\ell]$ is discretized at the translated extrema nodes of Chebychev
polynomials defined by $x_{i}=\ell/2(1-\cos(i\pi/M))$ for $i=0,\dots,M$. Let
us define the state variables $\bm{\rho}^{(m)}_{i}(t)=\rho^{(m)}(t,x_{i})$ and
$\bm{\varphi}_{i}(t)=\varphi(t,x_{i})$, along with the state vectors
$\bm{\rho}^{(m)}=(\bm{\rho}^{(m)}_{0},\dots,\bm{\rho}^{(m)}_{M})^{T}$ and
$\bm{\varphi}=(\bm{\varphi}_{0},\dots,\bm{\varphi}_{M})^{T}$. It can be shown
(e.g., see Ascher and Greif (2011)) that the spatial derivative of
$\rho^{(m)}(t,x)$ evaluated at $x=x_{i}$ is approximately equal to the
$(i+1)$-th column of $\bm{D}\bm{\rho}^{(m)}(t)$, where $\bm{D}$ is the
$(M+1)\times(M+1)$ differentiation matrix defined component-wise by

$\displaystyle\bm{D}_{ij}$ $\displaystyle=$
$\displaystyle\sum_{n\not=j}\frac{1}{x_{j}-x_{n}},\qquad(i=j),$
$\displaystyle\bm{D}_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{x_{j}-x_{i}}\prod_{n\neq
i,j}\frac{x_{i}-x_{n}}{x_{j}-x_{n}},\qquad(i\not=j).$

Similar approximations are made for spatial derivatives of other flow
variables. By discretizing the PDEs in (1)-(2) at the collocation nodes, we
obtain the system of ODEs in the vector variables $\bm{\rho}^{(1)}$,
$\bm{\rho}^{(2)}$, and $\bm{\varphi}$ given by

$\displaystyle\frac{d\bm{\rho}^{(m)}}{dt}+\bm{D}\left(\frac{\bm{\rho}^{(m)}}{\bm{\rho}^{(1)}+\bm{\rho}^{(2)}}\odot\bm{\varphi}\right)$
$\displaystyle=$ $\displaystyle 0,$ (6)
$\displaystyle\bm{D}\left(\sigma^{2}_{1}\bm{\rho}^{(1)}+\sigma^{2}_{2}\bm{\rho}^{(2)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda}{2D}\frac{\bm{\varphi}\odot|\bm{\varphi}|}{\bm{\rho}^{(1)}+\bm{\rho}^{(2)}},\quad$
(7)

where $\odot$ denotes the component-wise product, and quotients are defined
component-wise as well. The boundary conditions are integrated into the ODEs
with the substitutions $\bm{\rho}^{(m)}_{0}(t)=s^{(m)}(t)$ and
$\bm{\varphi}_{M}(t)=\overline{q}$. The solutions of (6)-(7) are obtained
numerically with Matlab using the stiff ODE solver `ode15s` (Shampine and
Reichelt (1997)) using the collocation points $\\{x_{i}\\}$ corresponding to
$M=16$. A nontrivial mass matrix is included in the obvious way.

##  3 Method to Quantify Chaos in Pipe Flow

As discussed above, the periodic forcing function $\gamma(t)$ that determines
partial pressures at the pipe inlet according to (5) is the only time-varying
parameter. We decompose $\gamma$ into its frequency components and analyze the
effect that each individual component has on the solution. Let us define a
sinusoidal variation in hydrogen mass fraction as

$\gamma(t)=\overline{\gamma}\left(1+\kappa\sin(2\pi\omega t)\right),$ (8)

where $\omega$ is the frequency of the sinusoid, $\kappa$ is the amplitude
factor, and $\overline{\gamma}$ is the mean concentration of hydrogen about
which the sinusoid oscillates. Our previous study (Baker et al. (2023b))
examined partitions of the boundary parameter plane $(\omega,\kappa)$ that
resulted in monotonic, periodic, and chaotic response regions. Below, we
design a controller to suppress chaotic responses and eliminate the area in
the $(\omega,\kappa)$ plane that defines the chaotic response region. We
suppose that a control design is successful at suppressing chaos if it results
in an empty chaotic response region when applied. Before we present the
control formula in the following section, we first discuss how chaos in
solutions to the reduced system (6)-(7) can be rigorously quantified.

We consider a solution to be chaotic if it is highly sensitive to initial
conditions (Lorenz (1995)). The extent of chaos in a finite-dimensional system
can be quantified by the largest Lyapunov exponent of the system (Benettin et
al. (1980)). This measure provides an estimate on the exponential rate of
divergence between two solutions that begin their evolution relatively close
to one another. Numerical methods have been developed to approximate the
largest Lyapunov exponent of a finite-dimensional system (Wolf et al. (1985);
Rosenstein et al. (1993); Brown et al. (1991)). These methods require an
appropriate embedding dimension, which could be problematic for solutions that
are generated by an underlying infinite-dimensional PDE system. Therefore, we
measure the extent of chaos in the reduced order system by using a rate of
divergence between two specific solutions. Note that this measure does not
necessarily provide an estimate of the largest Lyapunov exponent of the
continuum system. However, because the discretized system (6)-(7) is an
approximation of the continuum system (1)-(2) that retains the strongly
dissipative term on the right-hand side, we suppose that observations of chaos
in the former system indicate chaos in the latter.

Consider two numerical solutions $\psi_{1}(t_{n})=\psi(t_{n},\ell)$ and
$\psi_{2}(t_{n})=\psi(t_{n},\ell)$ with $|\psi_{2}(0)-\psi_{1}(0)|<\delta$,
where $\delta>0$ is small relative to $\psi_{1}(0)$. Here, $\psi(t_{n},\ell)$
represents the value of any one of the numerical flow variables evaluated at
the pipe outlet $x=\ell$ at time $t_{n}=(n/N)T$ for $n=0,\dots,N$. As in our
previous study (Baker et al. (2023b)), the exponential divergence between the
two solutions is quantified by

$\displaystyle\mathcal{C}_{\psi}$ $\displaystyle=$
$\displaystyle\frac{1}{|I_{T}|}\sum_{t_{n}\in
I_{T}}\log|\psi_{2}(t_{n})-\psi_{1}(t_{n})|$ (9)
$\displaystyle\qquad-\frac{1}{|I_{0}|}\sum_{t_{n}\in
I_{0}}\log|\psi_{2}(t_{n})-\psi_{1}(t_{n})|,\quad$

with $I_{0}=[t_{n_{0}},t_{n_{1}}]$ and $I_{T}=[t_{n_{2}},t_{n_{3}}]$, where
$n_{0}<n_{1}<n_{2}<n_{3}$. If one flow variable shows chaotic behavior, then
we would expect the other flow variables to show similar chaotic behavior.
Because the measures $\mathcal{C}_{\rho^{(1)}}$, $\mathcal{C}_{\rho^{(2)}}$,
and $\mathcal{C}_{p}$ are each sensitive to parameters and scalings, we define
the _chaos measure_

$\displaystyle\mathcal{C}=\min(\mathcal{C}_{\rho^{(1)}},\mathcal{C}_{\rho^{(2)}},\mathcal{C}_{p}).$
(10)

Because of the sensitivity of individual measures, we suppose that the minimum
of these measures provides more sensitivity to the presence of chaos than
using any one individually. Large and positive $\mathcal{C}$ indicates
exponential divergence between the two solutions over some time interval.
However, small and positive $\mathcal{C}$ does not necessarily imply that the
solutions diverge from one another. The _chaotic interface_ in the plane
$(\omega,\kappa)$ is defined to be the function $\kappa_{*}(\omega)$, where
for each $\omega$, $\kappa_{*}(\omega)$ is the lower bound that satisfies
$\mathcal{C}(\omega,\kappa)>C$ for all $\kappa>\kappa_{*}(\omega)$. We find
that $C=0.5$ is a useful threshold to identify chaos for the system parameters
considered in this paper, and we use this threshold moving forward. We say
that the solution corresponding to the forcing parameter pair
$(\omega,\kappa)$ is chaotic if $\kappa>\kappa_{*}(\omega)$.

##  4 Control Design to Suppress Chaos

Several different control designs varying in complexity may each demonstrate a
capability of suppressing chaos. The design that we propose is based on
proportional feedback. Our motivation for feedback control comes from
industrial accessibility. That is, because pipelines have sensors that record
time series measurements of pressure and mass flux values at high frequency
rates, we suppose that a chaos suppression controller could be integrated into
a pipeline control system. We propose a control action to stabilize mass flow
defined by

$u(t)=\overline{\mu}-\mu\left(\varphi(t,0)-\varphi(0,0)\right),$ (11)

where $\overline{\mu}$ and $\mu$ are scalars that represent the baseline
compressor ratio and the gain of the proportional feedback controller,
respectively. When $\mu=0$, the feedback controller is deactivated and the
resulting solution is called the baseline solution. We say that the controller
with gain $\mu$ suppresses chaos at the operating point $(\omega,\kappa)$ if
$\mathcal{C}(\omega,\kappa;\mu)<0.5$, where $\mathcal{C}$ is defined as in the
case without feedback control ($\mu\equiv 0$). The chaotic interface,
$\kappa_{*}(\omega;\mu)$, corresponding to the gain $\mu$ is also defined as
in the case without feedback control using equations (9) and (10). Although
the control design in (11) is conceptually simple, we will demonstrate that it
is highly effective at suppressing chaotic dynamics over the entire
$(\omega,\kappa)$ parameter plane.

Figure 1: Suppressing chaos by increasing the control gain $\mu$ in (11). Flow
variables in phase space are evaluated at the pipe outlet $x=\ell$ for
$t\in[75,100]$. The concentration boundary condition parameters are
$\omega=0.5$ cyc/hr, $\kappa=0.85$, and $\overline{\gamma}=0.2$. The remaining
parameters are $\overline{\mu}=1.75$, $\overline{p}=4$ MPa, $\overline{q}=75$
kg/m2s, $T=100$ hr, $\ell=50$ km, $D=0.5$ m, $\lambda=0.011$, $\sigma_{1}=338$
m/s, and $\sigma_{2}=4\sigma_{1}$.

Figure 2: Suppressing period doubling bifurcations. The concentration
parameters are $\omega=0.1$ cyc/hr, $\kappa=0.98$, and
$\overline{\gamma}=0.2$. The remaining parameters are the same as those in
Fig. 1.

Let us consider an example. We select a concentration variation parameter pair
$(\omega,\kappa)=(0.5,0.85)$ that results in a chaotic solution when the
controller is deactivated, i.e. $\mathcal{C}(\omega,\kappa;0)>0.5$. In Fig. 1,
this chaotic baseline solution along with two additional simulations that
correspond to two nonzero control gains are displayed to demonstrate a
suppression of chaos as the control gain increases. All of the pipeline and
boundary condition parameters, except for the control gains, are unchanged
between the three simulations. The phase space portrait of the baseline
solution corresponding to $\mu=0$ is shown on the left-hand side of Fig. 1.
Uninfluenced by the controller, the solution in this case regularly changes
direction in phase space with sharp gradients. This type of behavior is
typical for uncontrolled solutions (Baker et al. (2023b)). In the middle of
Fig. 1, we show the control action for a gain of $\mu=0.0025$ and the
resulting phase portrait. By inspection, we find that the controlled solution
for this gain value exhibits less hysteresis than the baseline (uncontrolled)
solution, but still does not have a closed periodic orbit in the physical
phase space. Increasing the gain further to $\mu=0.006$, we obtain a periodic
phase portrait with smooth gradients and a stabilizing periodic control
associated with this solution, as shown on the right-hand side of Fig. 1.
While this example demonstrates a transition to periodicity as the gain
increases, this does not show quantification of chaos in the rigorous
mathematical sense, i.e. that $\mathcal{C}(\omega,\kappa;\mu)<0.5$ for one of
the nonzero values of $\mu$. Nevertheless, one of the results in our previous
study (Baker et al. (2023b)) shows that the areas in the $(\omega,\kappa)$
parameter plane that correspond to periodic and non-chaotic dynamics in fact
coincide. Therefore, we associate the transition to periodicity with the
suppression of chaos in a manner that is sufficient to confirm by visual
inspection of simulations.

Although the above example may demonstrate a capability of the controller to
suppress chaos and stabilize downstream flows to a periodic orbit for the
operating point $(\omega,\kappa)=(0.5,0.85)$, the effectiveness of the
controller does not appear to be limited to a particular operating condition.
For another example, consider the operating frequency $\omega=0.1$ cyc/hr. The
baseline solution associated with this frequency approaches a periodic orbit,
as seen on the left-hand side of Fig. 2; however, this baseline solution
experiences a sequence of period doubling bifurcations as the amplitude factor
$\kappa$ increases to unity. For the amplitude factor $\kappa=0.98$, the
baseline solution has twice the period of the forcing concentration
$\gamma(t)$ as a result of the two loops that comprise the orbit, as seen on
the left-hand side of Fig. 2. Inspecting Fig. 2, we observe that the
controller suppresses period doubling, in the sense that the two loops in the
orbit of the baseline solution transition to a single loop for a sufficiently
large control gain. The examples in Figs. 1-2 demonstrate that the controller
may have the ability to suppress chaos while simultaneously stabilizing its
dynamics to a simple periodic orbit with the use of sufficiently large control
gains that depend on the operating points (and other pipeline parameters).
Moreover, the examples illustrate that a stabilizing controller has the
potential to significantly reduce large and undesirable gradients in phase
space that could otherwise appear in the uncontrolled solution.

Figure 3: Chaotic interfaces $\kappa_{*}(\omega;\mu)$ in the $(\omega,\kappa)$
plane as functions of $\omega$ for various control gains $\mu$. The region
defined by $0\leq\omega\leq 2$ and $0.5\leq\kappa\leq 1$ is discretized into a
$21\times 15$ grid of pairs of points. For each pair and for each control
gain, two solutions are simulated for 100 hours whose initial conditions
correspond to withdrawal rates $\overline{q}=75$ and $\overline{q}=75.1$. The
value $\mathcal{C}(\omega,\kappa;\mu)$ is determined from the two solutions
using (10) with intervals $I_{0}=[0.08N,0.15N]$ and $I_{T}=[0.5N,0.8N]$, where
$N=10,000$. The chaotic interface curve $\kappa_{*}(\omega;\mu)$ is obtained
from the set of values $\mathcal{C}(\omega,\kappa;\mu)$ on the grid according
to its definition and then interpolated with shape-preserving splines. The
remaining parameters are the same as those in Fig. 1.

We conclude by analyzing the chaotic interface $\kappa_{*}(\omega;\mu)$ as a
function of frequency for various control gains. Figure 3 depicts three
chaotic interfaces in the $(\omega,\kappa)$ parameter space that correspond to
the three control gains listed in the figure. It may appear from Fig. 3 that
for any two chaotic interfaces, $\kappa_{*}(\omega;\mu_{1})$ and
$\kappa_{*}(\omega;\mu_{2})$, the inequality
$\kappa_{*}(\omega;\mu_{1})\leq\kappa_{*}(\omega;\mu_{2})$ is satisfied for
all $\omega$ whenever $\mu_{1}\leq\mu_{2}$. However, from extensive numerical
simulations, we note that there are cases for which this proposed inequality
does not hold. Such non-monotonicity in computed criteria may be a consequence
of sensitivity issues that arise in characterizing chaotic dynamics
numerically. Nonetheless, the chaotic interfaces shown in Fig. 3 are in
agreement with the conclusions from the examples presented above. The
interfaces successfully demonstrate that a sufficiently large control gain
$\mu$ in a proportional feedback controller of the form in (11) could in fact
be capable of suppressing chaos over the entire ensemble of operating points
that comprise the boundary parameter plane.

##  5 Conclusions

We have reviewed a new problem in the boundary control of chaotic dynamics in
a distributed-parameter system defined by partial differential equations. This
setting has a compelling application to suppressing the chaotic flows that can
potentially occur in the downstream pressure of a natural gas pipeline into
which a time-varying mass fraction of hydrogen gas is injected. We demonstrate
a simple single-pipe test system in which the undesirable chaotic dynamics
emerge, and then propose a local controller that modulates the pipe inlet
pressure in order to stabilize the inlet flow using proportional feedback.
This controller is shown to completely suppress chaotic dynamics in flows and
pressures throughout the pipe for periodic hydrogen mass fraction variations
in the entire phase plane of frequency and amplitude parameters. Existing
compressor station control systems already incorporate the capability to
stabilize through-flow rates, in addition to a variety of pressure control
settings. Applying the type of controller proposed in our study, natural gas
pipelines could accommodate injections of hydrogen produced using renewable
energy to facilitate decarbonization while using capital investments in
infrastructure for their planned lifetimes.

## References

  * Alipoor and Mazaheri (2016) Alipoor, A. and Mazaheri, K. (2016).  Combustion characteristics and flame bifurcation in repetitive extinction-ignition dynamics for premixed hydrogen-air combustion in a heated micro channel.  _Energy_ , 109, 650–663.
  * Ascher and Greif (2011) Ascher, U.M. and Greif, C. (2011).  _A first course on numerical methods_.  SIAM.
  * Baker et al. (2023a) Baker, L.S., Kazi, S.R., Platte, R.B., and Zlotnik, A. (2023a).  Optimal control of transient flows in pipeline networks with heterogeneous mixtures of hydrogen and natural gas.  In _2023 American Control Conference (ACC)_ , 1221–1228. IEEE.
  * Baker et al. (2023b) Baker, L.S., Kazi, S.R., and Zlotnik, A. (2023b).  Transitions from monotonicity to chaos in gas mixture dynamics in pipeline networks.  _PRX Energy_ , 2, 033008.
  * Benettin et al. (1980) Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J.M. (1980).  Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory.  _Meccanica_ , 15, 9–20.
  * Boccaletti et al. (2000) Boccaletti, S., Grebogi, C., Lai, Y.C., Mancini, H., and Maza, D. (2000).  The control of chaos: theory and applications.  _Physics reports_ , 329(3), 103–197.
  * Botros and Jensen (2022) Botros, K.K. and Jensen, L. (2022).  Performance of twelve different equations of state for natural gas and hydrogen blends.  In _International Pipeline Conference_ , volume 86564, V001T08A002. American Society of Mechanical Engineers.
  * Brown et al. (1991) Brown, R., Bryant, P., and Abarbanel, H.D.I. (1991).  Computing the lyapunov spectrum of a dynamical system from an observed time series.  _Physical review A_ , 43(6), 2787.
  * Chaczykowski et al. (2018) Chaczykowski, M., Sund, F., Zarodkiewicz, P., and Hope, S.M. (2018).  Gas composition tracking in transient pipeline flow.  _Journal of Natural Gas Science and Engineering_ , 55, 321–330.
  * Dudnik et al. (1983) Dudnik, E.N., Kuznetsov, Y.I., Minakova, I., and Romanovskii, Y.M. (1983).  Synchronization in systems with strange attractors.  _Moscow Univ. Phys. Bull. Ser_ , 3(24), 84–87.
  * Eames et al. (2022) Eames, I., Austin, M., and Wojcik, A. (2022).  Injection of gaseous hydrogen into a natural gas pipeline.  _International Journal of Hydrogen Energy_ , 47(61), 25745–25754.
  * Eisenman et al. (2005) Eisenman, I., Yu, L., and Tziperman, E. (2005).  Westerly wind bursts: Enso’s tail rather than the dog?  _Journal of Climate_ , 18(24), 5224–5238.
  * Elaoud et al. (2017) Elaoud, S., Hafsi, Z., and Hadj-Taieb, L. (2017).  Numerical modelling of hydrogen-natural gas mixtures flows in looped networks.  _Journal of Petroleum Science and Engineering_ , 159, 532–541.
  * Götz et al. (2016) Götz, M., Lefebvre, J., Mörs, F., Koch, A.M., Graf, F., Bajohr, S., Reimert, R., and Kolb, T. (2016).  Renewable power-to-gas: A technological and economic review.  _Renewable energy_ , 85, 1371–1390.
  * Hafsi et al. (2019) Hafsi, Z., Elaoud, S., and Mishra, M. (2019).  A computational modelling of natural gas flow in looped network: Effect of upstream hydrogen injection on the structural integrity of gas pipelines.  _Journal of Natural Gas Science and Engineering_ , 64, 107–117.
  * Hu et al. (1995) Hu, G., Qu, Z., and He, K. (1995).  Feedback control of chaos in spatiotemporal systems.  _International Journal of Bifurcation and Chaos_ , 5(04), 901–936.
  * Kocarev (2001) Kocarev, L. (2001).  Chaos-based cryptography: a brief overview.  _IEEE Circuits and Systems Magazine_ , 1(3), 6–21.
  * Lorenz (1995) Lorenz, E.N. (1995).  _The essence of chaos_.  University of Washington press.
  * Melaina et al. (2013) Melaina, M.W., Antonia, O., and Penev, M. (2013).  Blending hydrogen into natural gas pipeline networks: A review of key issues.  Technical report.
  * Ott et al. (1990) Ott, E., Grebogi, C., and Yorke, J.A. (1990).  Controlling chaos.  _Physical review letters_ , 64(11), 1196.
  * Ozturk and Dincer (2021) Ozturk, M. and Dincer, I. (2021).  A comprehensive review on power-to-gas with hydrogen options for cleaner applications.  _International Journal of Hydrogen Energy_ , 46(62), 31511–31522.
  * Pizza et al. (2008) Pizza, G., Frouzakis, C.E., Mantzaras, J., Tomboulides, A.G., and Boulouchos, K. (2008).  Dynamics of premixed hydrogen/air flames in microchannels.  _Combustion and Flame_ , 152(3), 433–450.
  * Pyragas (1992) Pyragas, K. (1992).  Continuous control of chaos by self-controlling feedback.  _Physics letters A_ , 170(6), 421–428.
  * Rosenstein et al. (1993) Rosenstein, M.T., Collins, J.J., and De Luca, C.J. (1993).  A practical method for calculating largest lyapunov exponents from small data sets.  _Physica D: Nonlinear Phenomena_ , 65(1-2), 117–134.
  * Schuster and Stemmler (1997) Schuster, H.G. and Stemmler, M.B. (1997).  Control of chaos by oscillating feedback.  _Physical Review E_ , 56(6), 6410.
  * Shampine and Reichelt (1997) Shampine, L.F. and Reichelt, M.W. (1997).  The matlab ode suite.  _SIAM journal on scientific computing_ , 18(1), 1–22.
  * Subani et al. (2017) Subani, N., Amin, N., and Agaie, B.G. (2017).  Leak detection of non-isothermal transient flow of hydrogen-natural gas mixture.  _Journal of Loss Prevention in the Process Industries_ , 48, 244–253.
  * Tziperman et al. (1994) Tziperman, E., Stone, L., Cane, M.A., and Jarosh, H. (1994).  El Niño chaos: Overlapping of resonances between the seasonal cycle and the pacific ocean-atmosphere oscillator.  _Science_ , 264(5155), 72–74.
  * Wolf et al. (1985) Wolf, A., Swift, J.B., Swinney, H.L., and Vastano, J.A. (1985).  Determining lyapunov exponents from a time series.  _Physica D: nonlinear phenomena_ , 16(3), 285–317.
  * Zhang et al. (2009) Zhang, H., Liu, D., and Wang, Z. (2009).  _Controlling chaos: suppression, synchronization and chaotification_.  Springer Science & Business Media.
  * Zlotnik et al. (2023) Zlotnik, A., Kazi, S.R., Sundar, K., Gyrya, V., Baker, L., Sodwatana, M., and Brodskyi, Y. (2023).  Effects of hydrogen blending on natural gas pipeline transients, capacity, and economics.  In _PSIG Annual Meeting_ , PSIG–2312. PSIG.