On the Beta Transformation
Linas Vepstas







On the Beta Transformation
    Linas Vepstas
    December 2017 (Updated Feb 2018 and Dec 2018)
=================================================



linasvepstas@gmail.com doi:10.13140/RG.2.2.17132.26248

The beta transformation is the iterated map β x 1. The
special case of β=2 is known as the Bernoulli map, and is exactly
solvable. The Bernoulli map provides a model for pure, unrestrained
chaotic (ergodic) behavior: it is the full invariant shift on the
Cantor space { 0,1} ^ω. The Cantor space consists
of infinite strings of binary digits; it is notable for many properties,
including that it can represent the real number line. 

The beta transformation defines a subshift: iterated on the unit interval,
it singles out a subspace of the Cantor space that is invariant under
the action of the left-shift operator. That is, lopping off one bit
at a time gives back the same subspace.

The beta transform seems to capture something basic about the multiplication
of two real numbers: β and x. It offers insight into the
nature of multiplication. Iterating on multiplication, one would get
β^nx – that is, exponentiation; the mod 1 of the beta transform
contorts this in strange ways.

Analyzing the beta transform is difficult. The work presented here
is more-or-less a research diary: a pastiche of observations and some
shallow insights. One is that chaos seems to be rooted in how the
carry bit behaves during multiplication. Another is that one can surgically
insert “islands of stability” into chaotic (ergodic) systems,
and have some fair amount of control over how those islands of stability
behave. One can have islands with, or without a period-doubling “route
to chaos”.

The eigenvalues of the transfer operator seem to lie on a circle of
radius 1/β in the complex plane. Given that the transfer operator
is purely real, the appearance of such a quasi-unitary spectrum unexpected.
The spectrum appears to be the limit of a dense set of quasi-cyclotomic
polynomials, the positive real roots of which include the Golden and
silver ratios, the Pisot numbers, the n-bonacci (tribonacci, tetranacci,
etc.) numbers. 




§ INTRODUCTION


The last three or four decades of mathematical research has seen dramatic
advances in the theory of subshifts. This text is mostly not about
that, except to point out that this theory has very broad and important
impact on many branches of physics and mathematics. From the perspective
of the amateur enthusiast, the theory of subshifts finally exposes
and makes clear some of the mysterious and intriguing behavior of
fractals and of chaotic dynamical systems.

This text focuses almost entirely on just one particular map of the
unit interval, the β-transform, defined as the iterated map
β x1. As such, it is an example of an iterated map on the
unit interval of the real number line. Such maps have the form 

    f:[0,1]→[0,1]

and the topic is the exploration of the consequence of iterating the
map by composing: 

    f^n(x)=(f∘ f∘⋯∘ f)(x)=f(f(⋯ f(x)⋯))

Such one-dimensional iterated maps have been heavily studied, and
there is a large body of results, interconnecting many different concepts
and results from mathematics, and so having a particularly broad range. 

This text attempts to report some brand-new results on the β-transform.
This is perhaps surprising, as one might think that the β-transform
is sufficiently simple so as to be well-studied and well-understood,
it being among the very simplest of iterated one-dimensional maps.
This text also attempts to report these results in a naive and unsophisticated
fashion, in the hope that this makes the text readable for the interested
student and casual enthusiast.

Thus, although the author is personally excited by the advances in
the field, this text is neither a survey of known results on the β-transform,
nor does it much glance at most of the typical avenues that are available
for studying one-dimensional maps. This text does focus extensively
on the spectrum of the transfer operator (the “Ruelle Perron Frobenius
operator”), and thus it contributes to the general “Koopmania”.
Little prior knowledge is assumed, and the needed concepts are introduced
in a very casual and informal way. This will, no doubt, completely
discourage and dismay the formally trained mathematician. The best
I can offer is to reiterate: “new results”, off the beaten track.

This text begins with some pretty pictures, showing the iterated tent
and logistic maps, so as to remind the reader as to why this is an
interesting problem to study. The fact is that the β-transformation
is far more dry and abstract than the rather sexy logistic map, or
its complex cousin, the Mandelbrot set. The hope is that the β-transformation
is also simpler, and therefore, perhaps, easier to understand. The
reader will soon discover that there is nothing particularly easy
about it, and that, at every turn, one bumps into other interesting
areas of mathematics that could, perhaps should shed some light, but
don't actually seem to do so, in practice.

The most fun for the casual reader might be chapter 5, on the periodic
orbits, where the quasi-cyclotomic polynomials appear; these are polynomials
of the form p_n(z)=z^k+1-b_0z^k-b_1z^k-1-⋯-b_k-1z-1
for the b_j being binary bits (zero or one). Also quite fun is
the section on the islands of stability, which sheds light on how
one can take a purely ergodic (chaotic) system, and surgically insert,
as desired, islands of stability. The point here is that the classic
logistic map attracted interest precisely because of its interleaving
of chaos and stability; it turns out, one can manufacture such systems,
at will.

Perhaps the most surprising aspect of iteration (surprising to me)
is that the invariant measure consists of flat plateaus. Or perhaps
this is not surprising; after all, for the full shift, β=2,
the invariant measure is a straight line, and this is equivalent to
saying that the real numbers are aboslutely evenly distributed, or
that the toss of a perfectly fair coin is a perfectly even, random
sequences of heads and tails. Flatness implies an even distribution.
So flatness is not surprising. But look at it a different way: the
digit sequences are wildly and ergodically different, yet none-the-less
manage to sum to exactly the same value, when summed in an
analytic series! Why should crazy-different, crazy-chaotic sequences
yeild identical summations? That they do should be considered to be
the “fundamental theorem of analytic ergodics”. (I use the word
“analytic” here in the same sense as “analytic number theory”
or “analytic combinatorics”. Given a sequence of values, one explores
the analytic properties of its generating functions.) If this theorem
has ever been clearly stated before, or if it has been proven (for
any case other than the Bernoulli shift), I do not know. It seems
terribly important for analyzing subshifts.

A word about the format of this paper: this is a de facto “research
diary”, not a formal report. This, it contains various unfinished
corners and notes-to-self. 



 §.§ Bernoulli shift


The Bernoulli shift (aka the bit-shift map) is an iterated map on
the unit interval, given by iteration of the function

    b(x)=
    2x    0≤ x<1/2
    
    2x-1    1/2≤ x≤1

The symbolic dynamics of this map gives the binary digit expansion
of x. That is, write 

    b^n(x)=(b∘ b∘⋯∘ b)(x)=b(b(⋯ b(x)⋯))

to denote the n-fold iteration of b and let b^0(x)=x. The
symbolic dynamics is given by the bit-sequence

    b_n(x)=
    0    0≤ b^n(x)<1/2
    
    1    1/2≤ b^n(x)≤1

Of course, the symbolic dynamics recreates the initial real number:

    x=∑_n=0^∞b_n(x)2^-n-1

All of this is just a fancy way of saying that a real number can be
written in terms of it's base-2 binary expansion. That is, the binary
digits for x are the b_n=b_n(x), so that

    x=0.b_0b_1b_2⋯


The Bernoulli shift has many interesting properties, connecting it
to the Cantor set and to many self-similar fractals. I have explored
these in many other texts, as have other authors, and will not repeat
these here. The author is too lazy to provide a bibliography; the
reader is directed at search engines. 

The current task is to attempt to see how many of these properties
still hold in slightly more complex systems, and whether any of the
tools used to analyze and solve the Bernoulli shift can be applied
to these systems.



 §.§ Shift space


The use of the word “shift” here deserves a small bit of formality.
A “shift space” can be formally defined to be a set of infinite
sequences of a set of N letters (or symbols), together with a shift
operator T that takes each sequence, and lops off the left-most
symbol. For the Bernoulli shift, there are N=2 letters, taken from
the set { 0,1}. For the Bernoulli shift, one is typically
interested in the set of all possible infinite sequences: this is
the “full shift”. One writes { 0,1} ^ω
for this shift space, ω denoting countable infinity. For the
Bernoulli shift, the map b(x) is the shift operator:
it just lops of the left-most symbol.

In general, a shift space does not have to include every possible
sequence of symbols; it does, however, by definition, have to be shift-invariant.
That is, given some set S of infinite sequences of N symbols,
the set S is a shift space if and only if, by loping off the leading
symbol of each string, one regains S again. In formulas, a shift
space S must obey 

    TS=S


For example, S={ 000⋯,111⋯} contains only
two elements: the string of all zeros, and the string of all ones;
loping off the leading digit just returns S again. In general,
shift spaces may contain a finite, or a countable, or an uncountable
number of elements. In general, one defines the “full shift” as
the space N^ω of all possible strings of N symbols. Subsets
that are shift spaces are called “subshifts”.

The words “symbolic dynamics” also deserve some mention: given
one specific sequence x out of the shift space, one can ponder
“where it goes to”, as one lops off a symbol at a time. This gives
the “symbolic dynamics” or the “point dynamics” of the sequence.
The “orbit” is defined as the set { T^mx | m≥0}–
that is, the set of all places that x goes to. There are several
possibilities: one is that x is a fixed point, so that Tx=x.
Another is that x is a repeating sequence of symbols, in which
case iteration repeats as well: T^mx=x holds whenever the repeat
length is m; this is a periodic orbit. Most importantly, there
usually uncountably many non-periodic sequences or orbits. That is,
the number of periodic orbits is always countable: one merely arranges
them in lexicographic order, and one is done. As Cantor famously demonstrated
(and Hilbert so carefully expanded on) this cannot be done for the
non-periodic orbits: they are uncountable. 

In what follows, the text will in general confine itself to uncountable
case. Periodic orbits exist, but will be ignored; in a certain strict
sense, they constitute a set of measure zero. A number of glosses
like this will be made: for example, the real numbers, and the Cantor
space { 0,1} ^ω are both uncountable; however,
they are not in one-to-one correspondence, as some real numbers can
have two different representations as bit sequences. Specifically,
these are the fractions (2n+1)/2^m for positive integers
m,n – they can be validly represented by bit-sequences ending in
all-zeros, or all-ones. There are countably many such fractions, termed
the dyadic fractions. For the most part, this difference between the
real number line, and the Cantor space will be ignored.



 §.§ Beta shift


The beta shift is similar to the Bernoulli shift, replacing the number
2 by a constant real-number value 1<β≤2. It can be defined
as 


    T_β(x)=β x    0≤ x<1/2
    β(x-1/2)    1/2≤ x≤1

This map, together with similar maps, is illustrated in figure <ref>
below. 

Just as the Bernoulli shift generates a sequence of digits, so does
the beta shift: write

    k_n=
    0    0≤ T_β^n(x)<1/2
    
    1    1/2≤ T_β^n(x)≤1

Given the symbolic dynamics, one can reconstruct the original value
whenever 1<β as

    x=k_0/2+1/β(k_1/2+1/β(k_2/2+1/β(k_3/2+1/β(⋯))))

That is, one clearly sees that T_β(x) acts as a shift on
this sequence:

    T_β(x)=k_1/2+1/β(k_2/2+1/β(k_3/2+1/β(k_4/2+1/β(⋯))))

In this sense, this shift is “exactly solvable”: the above provides
a closed-form solution for iterating and un-iterating the sequence.

Multiplying out the above sequence, one obtains the so-called “β-expansion”
of a real number x, namely the series

    x=1/2 ∑_n=0^∞k_n/β^n

That is, the bit-sequence that was extracted by iteration can be used
to reconstruct the original real number. Setting β=2 in eqn
<ref> gives the Bernoulli shift. Explicitly, one
has T_2(x)=b(x).

Unlike the Bernoulli shift, not every possible bit-sequence occurs
in this system. It is a subshift of the full shift: it is a subset
of { 0,1} ^ω that is invariant under the action
of T_β. This is explored in greater detail in a later section. 



 §.§ Associated polynomial


The iterated shift can also be written as a finite sum. This is noted
here; it will be useful in later sections. Observe that 

    T_β(x)=β(x-k_0/2)

and that

    T_β^2(x)=β^2x-β/2(β k_0+k_1)

and that 

    T_β^3(x)=β^3x-β/2(β^2k_0+β k_1+k_2)

The general form is then:

    T_β^p(x)=β^px-β/2∑_m=0^p-1k_mβ^p-m-1

Since the k_m depend on both β and on x, this is not
a “true” polynomial; however, it will be useful in analysis, later. 



 §.§ Density Visualizations


Why is the beta transform interesting to explore? This can be partly
illustrated with some graphs. Shown in figure <ref>
is the “bifurcation diagram” for the beta transform. It visualizes
the long-term dynamics of the beta shift. Comparing to the usual bifurcation
diagram, e.g. for the Feigenbaum logistic map (shown in figure <ref>)
one thing becomes immediately apparent: there are no actual “bifurcations”,
no “islands of stability”, no period-doubling regions. Although
there are periodic orbits, these form a set of measure zero: the iteration
produces purely chaotic motion for all values of β. Thus, the
beta transform provides a clean form of “pure chaos”,[Formal mathematics distinguishes between many different kinds of chaotic
number sequences: those that are ergodic, those that are weakly or
strongly Bernoulli, weakly or strongly mixing. The beta transform
is known to be ergodic,<cit.> weakly mixing<cit.>
and weakly Bernoulli.<cit.>] without the pesky “islands of stability” popping up intermittently.

The visualization of the long-term dynamics is done by generating
a histogram, and then taking the limit, as follows. One divides the
unit interval into a fixed sequence of equal-sized bins; say a total
of N bins, so that each is 1/N in width. Pick a starting x,
and then iterate: if, at the n'th iteration, one has that j/N≤ b_β^n(x)<(j+1)/N,
then increment the count for the j'th bin. After a total of M
iterations, let c(j;M) be the count in the j'th bin. This count
is the histogram. In the limit of a large number of iterations, as
well as small bin sizes, one obtains a distribution: 

    ρ(y;x)=lim_N→∞lim_M→∞c(j;M)/Mj/N≤ y<j+1/N

This distribution depends on the initial value x chosen for the
point to be iterated; a “nice” distribution results when one averages
over all starting points:

    ρ(y)=∫_0^1ρ(y;x) dx

Numerically, this integration can be achieved by randomly sampling
a large number of starting points. Observe that ρ(y) is a probability
distribution:

    1=∫_0^1ρ(x) dx

This probability distribution is an eigenstate of the transfer operator
for the beta transform; the definition of the transfer operator of
the beta transform is given later. Probability distributions are the
same thing as measures; this particular distribution is invariant
under iteration, and thus is often called the invariant measure, or
sometimes the Haar measure.

For each fixed β, one obtains a distinct distribution ρ_β(y).
The figure <ref> illustrates
some of these distributions. Note that, for β<1, the distribution
is given by ρ_β(y)=δ(y), a Dirac delta function,
located at y=0.



The general trend of the distributions, as a function of β,
can be visualized with a Feigenbaum-style “bifurcation diagram”,
shown in figure <ref>. This color-codes
each distribution ρ_β(y) and arranges them in a stack;
a horizontal slice through the diagram corresponds to ρ_β(y)
for a fixed value of β. The term “bifurcation diagram”
comes from its use to visualize the logistic map iterator.






 §.§ Tent Map


The tent map is a closely related iterated map, given by iteration
of the function

    v_β(x)=β x    0≤ x<1/2
    β(1-x)    1/2≤ x≤1

Its similar to the beta shift, except that the second arm is reflected
backwards, forming a tent. The bifurcation diagram is shown in figure
<ref>. Its is worth contemplating the similarities
between this, and the corresponding beta shift diagram. Clearly, there
are a number of shared features.






 §.§ Logistic Map


The logistic map is related to the tent map, and is given by iteration
of the function

    f_β(x)=2β x(1-x)

It essentially replaces the triangle forming the tent map with a parabola
of the same height. That is, the function is defined here so that
the the same value of β corresponds to the same height for
all three maps. Although the heights of the iterators have been aligned
so that they match, each exhibits rather dramatically different dynamics.
The β-transform has a single fixed point for β<1, and
then explodes into a fully chaotic regime above that. By contrast,
the logistic map maintains a single fixed point up to β=3/2,
where it famously starts a series of period-doubling bifurcations.
The onset of chaos is where the bifurcations come to a limit, at β=3.56995/2=1.784975.
Within this chaotic region are “islands of stability”, which do
not appear in either the β-transform, or in the tent map. The
tent map does show a period-doubling regime, but in this region, there
are no fixed points: rather, the motion is chaotic, but confined to
multiple arms. At any rate, the period doubling occurs at different
values of β than for the logistic map.

The bifurcation diagram is shown in figure <ref>.
Again, it is worth closely examining the similarities between this,
and the corresponding tent-map diagram, as well as the β-transform
diagram. Naively, it would seem that the general structure of the
chaotic regions are shared by all three maps. Thus, in order to understand
chaos in the logistic map, it is perhaps easier to study it in theβ-transform.



The general visual similarity between the figures <ref>,
<ref> and <ref>
should be apparent, and one can pick out and find visually similar
regions among these three illustrations. Formalizing this similarity
is a bit harder, but it can be done: there is a way to make all three
of these maps be “topologically conjugate” to one-another. This
is perhaps surprising to some readers, but is based on the observation
that the “islands of stability” in the logistic map are countable,
and are in one-to-one correspondence with certain “trouble points”
in the iterated beta transformation. These are in turn in one-to-one
correspondence with rational numbers. With a slight distortion of
the beta transformation, the “trouble points” can be mapped to
the islands of stability, in essentially the same way that “phase
locking regions” or “Arnold tongues” appear in the circle map.
But this is all for a later section, again, mentioned here only to
whet the appetite.



 §.§ Beta Transformation


After exactly one iteration of the beta shift , all initial points
β/2≤ x≤1 are swept up into the domain 0≤ x<β/2,
and never leave. Likewise, the range of the iterated beta shift is
0≤ x<β/2. Thus, an alternative representation of the beta
shift, filling the entire unit square, can be obtained by dividing
both x and y by β/2 to obtain the function

    t_β(u)=β u    0≤ u<1/β
    β u-1    1/β≤ u≤1

which can be written more compactly as t_β(x)=β x 1.
In this form, the function is known as the beta transform, and is
often called the β-transformation, presenting a typesetting
challenge to search engines when used in titles of papers. The orbit
of a point x in the beta shift is identical to the orbit of a point
u=2x/β in the beta transformation. Explicitly comparing to
the beta shift of eqn <ref>:

    T_β^n(x)=β/2t_β^n(2x/β)

The beta shift and the β-transformation are essentially “the
same function”; this text works almost exclusively with the beta
shift, and is thus idiosyncratic, as it flouts the much more common
convention of working with the β-transformation. There is no
particular technical reason for this; it is rather due to happenstance. 

After a single iteration of the tent map, a similar situation applies.
After one iteration, all initial points β/2≤ x≤1 are swept
up into the domain 0≤ x<β/2. After a finite number of iterations,
all points 0<x≤β(1-β/2) are swept up, so that
the remaining iteration takes place on the domain β(1-β/2)<x<β/2.
It is worth defining a “sidetent” function, which corresponds
to the that part of the tent map in which iteration is confined. It
is nothing more than a rescaling of the tent map, ignoring those parts
outside of the above domain that “wander away”. The sidetent is
given by

    s_β(u)=β(u-1)+2    0≤ u<β-1/β
    β(1-u)    β-1/β≤ u≤1

Performing a left-right flip on the side-tent brings it closer in
form to the beta-transformation. The flipped version, replacing u→1-u
is 

    f_β(u)=β u    0≤ u<1/β
    
    2-β u    1/β≤ u≤1


The tent map (and the flipped tent) exhibits fixed points (periodic
orbits; mode-locking) for the smaller values of β. These can
be eliminated by shifting part of the tent downwards, so that the
diagonal is never intersected. This suggests the “sidetarp”:

    a_β(u)=β u    0≤ u<1/β
    β(1-u)    1/β≤ u≤1

The six different maps under consideration here are depicted in figure
<ref>. It is interesting to compare
three of the bifurcation diagrams, side-by-side. These are shown in
figure <ref>.








 §.§ Beta Transformation Literature Review and References


The β-transformation, in the form of t_β(x)=β x 1
has been well-studied over the decades. The beta-expansion <ref>
was introduced by A. Renyi<cit.>, who demonstrates the existence
of the invariant measure. The ergodic properties of the transform
were proven by W. Parry<cit.>, who also shows that the system
is weakly mixing. 

An explicit expression for the invariant measure was independently
obtained by A.O. Gel'fond<cit.> and by W. Parry<cit.>,
as

    ν_β(y)=1/F∑_n=0^∞ε_n(y)/β^n

where ε_n is the digit sequence 

    ε_n(y)=
    0    t_β^n(1)≤ y
    
    1

and F is a normalization constant. 

In the same way that a dyadic rational p/2^n has two different
binary expansions, one ending in all-zeros, and a second ending in
all-ones, so one may also ask if and when a real number x might
have more than one β-expansion (for fixed β). In general,
it can; N. Sidorov shows that almost every number has a continuum
of such expansions!<cit.> 

Conversely, the “univoke numbers” are those values of β
for which there is only one, unique expansion for x=1. These are
studied by De Vries.<cit.>

The β-transformation has been shown to have the same ergodicity
properties as the Bernoulli shift.<cit.> The fact that the
beta shift, and its subshifts are all ergodic is established by Climenhaga
and Thompson.<cit.>

An alternative to the notion of ergodicity is the notion of universality:
a β-expansion is universal if, for any given finite string
of bits, that finite string occurs somewhere in the expansion. This
variant of universality was introduced by Erdös and Komornik<cit.>.
Its is shown by N. Sidorov that almost every β-expansion is
universal.<cit.> Conversely, there are some values of β
for which rational numbers have purely periodic β-expansions;<cit.>
all such numbers are Pisot numbers.<cit.>

The symbolic dynamics of the beta-transformation was analyzed by F.
Blanchard<cit.>. A characterization of the periodic points
are given by Bruno Maia<cit.>. A discussion of various open
problems with respect to the beta expansion is given by Akiyama.<cit.> 

When the beta expansion is expanded to the entire real-number line,
one effectively has a representation of reals in a non-integer base.
One may ask about arithmetic properties, such as the behavior of addition
and multiplication, in this base - for example, the sum or product
of two β-integers may have a fractional part! Bounds on the
lengths of these fractional parts, and related topics, are explored
by multiple authors.<cit.>

Certain values of β - generally, the Pisot numbers, generate
fractal tilings,<cit.>
which are generalizations of the Rauzy fractal. An overview, with
common terminology and definitions is provided by Akiyama.<cit.>
The tilings, sometimes called (generalized) Rauzy fractals, an be
thought of as living in a direct product of Euclidean and p-adic
spaces.<cit.>

The set of finite beta-expansions constitutes a language, in the formal
sense of model theory and computer science. This language is recursive
(that is, decidable by a Turing machine), if and only if β
is a computable real number.<cit.>

The zeta function, and a lap-counting function, are given by Lagarias<cit.>.
The Hausdorff dimension, the topological entropy and general notions
of topological pressure arising from conditional variational principles
is given by Daniel Thompson<cit.>. A proper background
on this topic is given by Barreira and Saussol<cit.>. 

None of the topics or results cited above are made use of, or further
expanded on, or even touched on below. That's mostly because what
follows below is just a ... different perspective on the topic. Not
intentionally so; it just happened to come out that way.



§ SYMBOLIC DYNAMICS


The Bernoulli shift corresponds to the sequence of binary digits of
a real number. Such sequences can be imagined to belong to the space
of all possible sequences of binary digits, the Cartesian product
of infinitely many copies of the set containing two elements:

    { 0,1}×{ 0,1}×{ 0,1}×⋯={ 0,1} ^ω=2^ω

This space has a natural topology, the product topology, which differs
sharply from the natural topology on the real-number line. Essentially
all of the strange phenomena of fractals and of iterated functions
follows from the product topology on this sequence. 

One notable effect that can be explained in terms of the product topology
is the fractal self-similarity of many kinds of fractals: this arises
from the self-similarity of the product space under the action of
a shift: specifically, the left-shift, which discards the left-most
digit, and shifts the rest over by one. The shift operator itself
is that operator that performs this shift; self-similar fractals can
be seen to be eigenstates of the shift operator. 

Another notable effect is the close proximity of the Cantor set to
the proceedings. In a certain sense, the Cantor set can be understood
to be the most basic manifestation of the product space. When attuned
to its presence, it can be seen everywhere throughout the proceedings.

A third byproduct is the manifestation of the infinitely-deep binary
tree. This arises when the set { 0,1} of the product
space is re-interpreted as the set { L,R} of left-right
moves. At each point in a binary sequence, one can make a choice of
one of two things: to move left or right. This naturally suggests
a binary decision tree. 

A fourth byproduct is the presence of some implicit, ambient hyperbolic
space. The infinite binary tree, when drawn on flat two-dimensional
space, simply “runs out of room”, as each subsequent branching
pushes closer together. The infinite binary tree can be embedded in
the simplest hyperbolic space, the Poincaré disk or upper-half-plane,
in such a way that the distance, the spacing between two neighboring
nodes is always the same. Visually, this takes the form of some prototypical
M.C. Escher drawing, of a repeated fractal form moving out to the
edge of a disk. This makes the self-similar shape of the infinite
binary tree manifest: as one moves from one location to another, one
always sees “the same thing” in all directions: the space is homogeneous.

The rational numbers play a very special role in the infinite binary
tree. Dyadic rationals, of the form (2p+1)/2^n for
integers p and n correspond to bit sequences (eqn <ref>)
that terminate in all-zeros after a finite number of moves. That is,
after an initial “chaotic” sequence, they settle down to a fixed
point of period one. General rational numbers p/q behave similarly,
in that after an initial “chaotic” sequence, they settle down
to periodic orbits of some fixed period. The bit-sequence becomes
cyclic. This cyclic behavior implies that most of classical number
theory can be dragged into the proceedings. Any particular statement
that classical number theory makes with regard to rational numbers,
or even modular forms, can be promptly ported over to a statement
about the bit-sequences and the orbits of the Bernoulli shift, usually
taking on a strange and unrecognizable form. 

All of these things go together, like hand in glove: whenever one
is manifest and visible, the others are lurking right nearby, in the
unseen directions. All of these things can be given a formal and precise
definition, and their explicit inter-relationships articulated. This
has been done by a wide variety of authors over the last four decades;
a proper bibliography would be overwhelming. I have written on all
of thee topics, trying to present them in the simplest, most jargon-free
way that I can, in a dozen different texts available wherever you
found this one. The ideas will not be repeated here; they are not
immediately useful to the current proceedings. None-the-less, the
general interplay between all of these concepts is extremely important
to understand, and burbles constantly under the surface of the current
proceedings. In essence, shifts and subshifts are interesting precisely
because they touch on so any different topics; and, conversely, so
many different areas of mathematics can inform the subshift. 



 §.§ Symbolic Dynamics


Given that iteration can generate strings of binary digits, and that
these can be reassembled back into real numbers, it is interesting
to ask what those mappings look like. To firm up the notation, let(b_n)=(b_0,b_1,,⋯)
denote a sequence of bits (or symbols) and write


    x_β((b_n))=1/2 ∑_n=0^∞b_n/β^n

as the real number generated from that sequence. Conversely, given
a real number x, let (k_n;β(x))
denote the sequence of bits obtained by iterating the beta shift on
x with constant β; that is, the sequence generated by eqn.
<ref>. The bit sequence for (k_n;2(x))
is just the bit sequence (b_n(x)) generated
by eqn <ref>. The transformations between symbol
sequences and real numbers make sense only when 1<β≤2.

Two interesting functions can be considered. One is the compressor

    _β(y)=x_2((k_n;β(y)))

and the other is the expander

    _β(y)=x_β((k_n;2(y)))

The terms “compressor” and “expander” are being invented here
to indicate negative and positive Lyapunov exponents associated with
the two functions. For almost all y, the compressor function is
pushing nearby points closer together; the total measure of the range
of the compressor function is less than one. Likewise, for almost
all y, the expander function is pushing nearby points apart. These
two functions are illustrated in figures <ref>
and <ref>.





The two functions are adjoint; specifically, one has that _β(_β(y))=y
but that _β(_β(y)) y.
The former relation is equivalent to eqn. <ref>.
Not all possible sequences of bit strings appear in the beta shift
sequence (k_n;β(x)); that is, this
function is not a surjection onto { 0,1} ^ω.
This manifests itself as the gaps in the range of the compressor function,
clearly visible in figure <ref>. If a sequence
of bits is viewed as a sequence of left-right moves walking down a
binary tree, this implies that some branches of the tree are never
taken, and can be pruned. Only branches on the right are ever pruned:
That is, there can be arbitrarily long sequences of zeros in the expansion,
but the longest possible sequence of 1's is always bounded. The longest
run of 1's possible is the largest value of n that satisfies 

    2≥1+β+β^2+⋯β^n-1/β^n-1

Solving, the bound is

    n=1+⌊-log(2-β)/logβ⌋

That is, every n'th right branch is pruned from the binary tree.
For example, a run of three 1's in a row is possible only for β≥(1+√(5))/2=1.618034⋯
the Golden Ratio. The range of _β(y)
is most of, but not all of the Cantor set. The figure <ref>
visualizes the range of the compressor as a function of β. 






 §.§ Shifts with holes


Viewed as a shift space, as opposed to a cut-down binary tree, the
trimming can be thought of as a punching of holes into the full shift.
This requires a bit of mental gymnastics. Let (a,c)
be an (open) interval on the real number line: (a,c)={ x | a<x<c}.
Given the Bernoulli shift b(x)=T_2(x) from
eqns <ref> or <ref>, consider the
set

    ℐ(a,c)={ x | b^n(x)∉(a,c)n≥0}

That is, as one iterates on some fixed x, one requests that no
iterate b^n(x) ever lands in the interval (a,c).
In essence, one has punched a hole in the unit interval; this corresponds
to a “hole” in the full Bernoulli shift. The set ℐ(a,c)
is what remains after punching such a hole. 

How can this be visualized? Considering the case n=0, its clear
that ℐ(a,c) cannot contain (a,c).
That is, ℐ(a,c)∩(a,c)=∅.
For n=1, the interval (a,c) can come from one of
two places: either from (a/2,c/2) or
from (a+1/2,c+1/2), and so neither of
these can be in ℐ(a,c). Continuing, for n=2,
the intervals (a/4,c/4), (a+1/4,c+1/4),
(a+2/4,c+2/4) and (a+3/4,c+3/4)
must also be gone. Continuing in this fashion, one proceeds with an
infinite hole-punch: to obtain ℐ(a,c), one
just cuts out (a,c) and everything that iterates to
(a,c). For the holes, write

    ℋ(a,c)=⋃_n=0^∞⋃_k=0^2^n-1(a+k/2^n,c+k/2^n)

and for the interval with the holes punched out: 

    ℐ(a,c)=[0,1]\ℋ(a,c)

where ⋃ denotes set-union and \ denotes set subtraction.
It is not hard to see that, in the end, this forms a contorted Cantor
set, using the standard midpoint-subtraction algorithm, but with different
endpoints. The canonical Cantor set is built by taking (a,c)=(1/3,2/3).

Note that both ℋ(a,c) and ℐ(a,c)
are subshifts: applying the left-shift to them just returns the same
set again. Bot are invariant under the action of the shift operator.
In formulas, 

    bℋ(a,c)=ℋ(a,c)

 and 

    bℐ(a,c)=ℐ(a,c)

where, for notational simplicity, the parenthesis are not written,
so that for the set S, write bS=b(S). As shifts,
its more appropriate to view both as sets of bit-sequences, so that
the proper relationship between one and the other should have been
written as

    ℐ(a,c)={ 0,1} ^ω\ℋ(a,c)


How should these subshifts be visualized as strings? Let (b_n(x))
be the bit sequence generated by x, for some a<x<c. The cut
operation states that such strings can never occur anywhere in ℐ(a,c).
Explicitly, ℐ(a,c) never contains sequences
of the form d_0d_1d_2⋯ d_kb_0(x)b_1(x)b_2(x)⋯
for any arbitrary leading bits d_0d_1d_2⋯ d_k. 

How should these subshifts be visualized as binary trees? The simplest
case to visualize is to take a=m/2^n and c=(m+1)/2^n
being dyadic rationals, for some integers m,n. In this case, one
takes the bit-expansion for both have the same n leading bits:
one starts at the root of the tree, and walks down the binary tree,
making left-right moves in accordance with this sequence, and after
n moves, arrives at a node above a subtree. Just cut out this subtree,
in it's entirety. That's the first cut. Now repeat the process, for
the left and right subtrees, from off the root, ad infinitum.
For a and c not dyadic rationals, the process is more complicated.
If a and c are ordinary rationals, thus having a repeating bit-sequence,
one performs in the same way, but cyclically walking down the side
branches of subtrees. For a and c irrational, the algorithm
is considerably more complicated, and is left as an exercise for the
reader :-).

A general classification of shifts with holes, for the beta transform,
was performed by Lyndsey Clark<cit.>.



 §.§ Generalized compressors and expanders


The range of the compressor function is a shift with a hole. Specifically,
for a given β, the range of _β is ℐ(β/2,1/2).
The construction for shifts with holes can then be applied to construct
generalized compressor and expander functions. One way, which is really
rather cheesy, but it works, is to define the function

    _β,γ(a;x)=∑_n=0^∞[1/γ^n+1∑_k=0^2^n-1δ(x-a+k/β^n)]

and then define the generalized compressor as

    (a;x)=∫_0^x(a;y)dy

That is, as one walks along the unit interval, from left to right,
one picks up points with weights on them, obtaining a generalized
Devil's staircase (Cantor-Vitali) function. This generalization does
not seem to be terribly useful here, and is left to rot.



 §.§ Self-similarity


Subshifts are, by definition, self-similar. If S is a subshift,
and T is the shift operator, then TS=S is a part of the definition
of the subshift. It is fun to see how this actually manifests itself
on the unit interval.

So, the two functions cpr and pdr are self-similar. The pdr function
demonstrates classic period doubling self-similarity: namely, under
g(x)=x/2, it behaves as

    (_β∘ g)(x)=_β(x/2)=1/β_β(x)

while under reflection r(x)=1-x, it behaves as 

    (_β∘ r)(x)=_β(1-x)=β/2(β-1)-_β(x)

Note that 

    lim_x→1 _β(x)=β/2(β-1)

The full dyadic monoid is generated by the generators g and r;
see other posts from me for lengthly expositions on the structure
of the dyadic monoid and its relationship to the Cantor set and a
large variety of fractals. 

Here, g is the generator that corresponds to the shift operator
T. The notation g is used only to stay consistent with other
things that I've written. The generator r indicates that the subshift
is also invariant under reflection; in this case, under the exchange
of the symbols 0↔1 in the corresponding shift.

The function cpr also exhibits self-similarity, although it alters
(expands) what happens on the x axis. Several self-similarities
are apparent. First, for 0≤ x≤1, one has

    _β(x/2)=1/2_β(β x/2)

Equivalently, for 0≤ y≤β/2 one can trivially restate the
above as

    _β(y/β)=1/2_β(y)

Although it follows trivially, this restatement helps avoid later
confusion.

The left and right halves are identical to one-another, but offset:

    _β(1/2+x/2)=1/2+_β(x/2)

It follows that


    _β(1/2+y/β)=1/2+1/2_β(y)

Combining the above results into one, one has that


    _β(y/β)+_β(1/2+y/β)=1/2+_β(y)

This last form is interesting, as it makes an appearance in relation
to the transfer operator, defined below.



 §.§ Other things with similar symmetry


The cpr curve is just one that belongs to a class of such curves.
As an example, one may construct a Takagi (blancmange) curve by iterating
triangles whose peak is located at 1/β. The Takagi curve is
an example of a curve transforming under a 3-dimensional representation
of the dyadic monoid; the cpr curves transforms under a two-dimensional
representation. See my paper on the Takagi curve for details. Figure
<ref> shows such a curve. Denote by _β;w(x)a
curve constructed in this fashion. The transformation properties of
this curve include self-similarity on the left, as

    _β;w(x/β)=x+w _β;w(x)

for 0≤ x≤1 and self-similarity on the right, as

    _β;w(1/β+x(1-1/β))=1-x+w _β;w(x)

Both of these properties follow directly from the construction of
the curve; they can be taken as the defining equations for the curve.
That is, the curve can be taken as that function which satisfies these
two recursion relations.

The derivative of the skew Takagi curve is shown in figure <ref>,
and, for lack of a better name, could be called the skew Haar fractal
wavelet. It can be defined as the formal derivative

    _β;w(x)=d/dx _β;w(x)

This formal derivative is well-defined, as the skew Takagi is smooth
and piecewise-linear almost everywhere; the places where it has corners
is a dense set of measure zero. That is, the derivative is defined
everywhere, except on a set of measure zero, which happens to be dense
in the unit interval.



Note that the Haar fractal wavelet is piece-wise constant everywhere.
It is constructed from a “mother wavelet” given by

    h_β(x)=β   0≤ x<1/β
    -β/β-1   1/β≤ x≤1

which is then iterated on to form the fractal curve _β;w(x).
The self symmetries are

    _β;w(x/β)=β+w _β;w(x)

 and

    _β;w(1/β+x(1-1/β))=-β/β-1+w _β;w(x)







 §.§ Fixed Points; Periodic Orbits


The Bernoulli shift, given by eqn <ref>, generates
every possible bit-sequence. As was observed in a previous section,
not every possible bit-sequence occurs in the beta shift. The longest
sequence of all-ones possible was given by eqn <ref>.
Arbitrary finite lengths of zeros do appear; but are there fixed points,
i.e. sequences that terminate in all-zeros? Clearly, x=1/2β^n
is such a fixed point: after n+1 iterations of eqn <ref>,
x goes to zero, and stays there. Is this the only such fixed point?
The answer depends on β. If β can be written in the
form of β^n=2m+1 for some integers n and m, then the
values of x which can iterate down to zero in n+1 steps are
dense in the interval [0,β/2]. Curiously, such values
β are dense in the interval [1,2). A later chapter
performs explores periodic orbits in great detail.



§ TRANSFER OPERATORS


The discovery and study of invariant measures, as well as of decaying
states can be approached via the transfer operator, or, properly named,
the Ruelle-Frobenius-Perron operator. This is an operator that captures
the behavior of a distribution under the action of a map. The invariant
measure is an eigenstate of this operator; indeed, it provides a formal
definition for what it means to be invariant under the action of the
map.

Given an iterated map g:[0,1]→[0,1] on the unit interval, the
transfer operator defines how distributions are acted on by this map.
It is defined as

    [ℒ_gf](y)=∑_x=g^-1(y)f(x)/|g^'(x)|

The left adjoint of the transfer operator is the composition operator
(Koopman operator). This is defined as 

    [𝒞_gf](y)=f(g(y))

The Koopman operator is adjoint, in the sense that ℒ_g𝒞_g=1
but that, in general, 𝒞_gℒ_g1. 



 §.§ The β-transform Transfer Operator


The transfer operator for the beta shift map T_β(x) is

    [ℒ_βf](y)=1/β[f(y/β)+f(y/β+1/2)]    0≤ y≤β/2
    
    0    β/2<y≤1

or, written more compactly

    [ℒ_βf](y)=1/β[f(y/β)+f(y/β+1/2)]Θ(β/2-y)

where Θ is the Heaviside step function. The density distributions
graphed in figure <ref> are those
functions satisfying 

    [ℒ_βρ](y)=ρ(y)

That is, the ρ(y) satisfies 

    ρ(y)=1/β[ρ(y/β)+ρ(y/β+1/2)]Θ(β/2-y)

This is generally referred to as the Ruelle-Frobenius-Perron (RFP)
eigenfunction, as it corresponds to the largest eigenvalue of the
transfer operator, and specifically, the eigenvalue 1. 

More generally, one is interested in characterizing the eigenspectrum

    [ℒ_βρ](y)=λρ(y)

for eigenvalues |λ|≤1 and eigenfunctions ρ(y).
Solving this equation requires choosing a space of functions in which
to work. Natural choices include any of the Banach spaces, and in
particular, the space of square-integrable functions. Particularly
interesting is the space of almost-smooth functions, those having
discontinuities at only a countable number of locations, but otherwise
being infinitely differentiable. Although the discussion so far implicitly
conditions one to restrict oneself to real-valued functions, and to
consider only real-valued eigenvalues, this is perhaps too sharp a
restriction. It will become clear in the following chapters that even
the most basic form of ℒ_β has a complex-valued
spectrum. At any rate, it should be obvious that, whatever the choice
of function space, one must have that ρ(y)=0 whenever
β<2y. This turns out to be a rather harsh condition.

At least one basic fact is known: for at least some kinds of function
spaces, the RFP eigenfunction is given by Gel'fond and Parry, as shown
in eqn <ref>. More precisely, it is just the
rescaled form ρ(x)=ν(2x/β) for x<β/2
and zero otherwise. Changing vocabulary, this is sometimes called
the “invariant measure”; as it describes a measure on the unit
interval. That is, for the space of all possible measures on the unit
interval, the Gel'fond-Parry measure is one of the eigenfunctions
of the transfer operator. Some caution is advised here: for the special
case of β=2, that is, the Bernoulli shift, one has as an eigenfunction
the Minkowski measure<cit.>; it has eigenvalue 1, but
is otherwise quite pathological: it is continuous nowhere, zero on
the rationals, and divergent on the rest (i.e. on a “fat” Cantor
set). There's no particular reason to think that this holds only for
β=2; measures can be, in general, very unusual functions.

A very minor simplification can be achieved with a change of variable.
Let y=β/2-ϵ. Then the eigenequation becomes

    λβρ(β/2-ϵ)=ρ(1/2-ϵ/β)+ρ(1-ϵ/β)

The second term vanishes whenever β/2<1-ϵ/β or ϵ<β(1-β/2)
and so one has the simpler recurrence relation

    λρ(y)=1/βρ(y/β)

whenever β(β-1)<2y≤β.

The equations <ref> and <ref> can be
treated as recurrence relations, defining the λ=1 eigenstate.
Recursing on these gives exactly the densities shown in figure <ref>.
Computationally, these are much, much cheaper to compute, at least
for β much larger than 1, although convergence issues present
themselves as β approaches 1. The resulting density may be
called the Ruelle-Frobenius-Perron eigenstate; because it can be quickly
computed, it provides an alternative view of figure <ref>,
free of stochastic sampling noise.



 §.§ A note on complex eigenvalues


Since the operator ℒ_β is purely real, then if
it has a complex spectrum, the eigenvalues and eigenfunctions must
come in complex-conjugate pairs. This can make numerical searches
and numerical convergence behave in unexpected ways, so some brief
commentary is in order. 

Assume that there exists some complex-valued eigenfunction ρ_λ(x)
for fixed, complex eigenvalue λ. Write it's real and complex
components as

    ρ_λ(x)=σ(x)+iτ(x)

while also writing λ=a+ib. Then

    [ℒ_βσ](x)=    1/2[ℒ_β(ρ_λ+ρ̅_λ)](x)
    
    =    1/2(λρ_λ(x)+λ̅ρ̅_λ(x))
    
    =     aσ(x)-bτ(x)

Both left and right sides of the above are real. If one had somehow
stumbled upon σ(x) numerically, as an eigenvector-candidate,
then the above admixing of the imaginary component would quickly throw
one off the hunt. Thus, a numeric search for complex-valued eigenfunctions
must necessarily take into account eigenfunction pairs, with real
and imaginary components that mix together as above.



 §.§ Almost-solutions


If one ignores the Heaviside step function in the definition <ref>,
one easily finds a number of “almost solutions” to the transfer
operator. These are most easily discussed by defining the operator


    [𝒫_βf](y)=1/β[f(y/β)+f(y/β+1/2)]

Solving this operator is relatively straight-forward. Consider, for
example, the monomial f(y)=y^n. Clearly, [𝒫_βy^n]
is a polynomial of degree n and that therefore, the space of polynomials
is closed under the action of 𝒫_β. But this result
is even stronger: the monomials provide a basis in which 𝒫_β
is upper-triangular, i.e. solvable. It's eigensolutions in
this basis are polynomials. The eigenspectrum is clearly discrete,
and is given by (β)^-n-1 for integers n corresponding
to the degree of the polynomial eigensolution. 

This all goes horribly wrong if one instead considers ℒ_β
and the almost-monomials f(y)=y^nΘ(β/2-y).
This does not provide a basis that is closed under the action of ℒ_β.
Attempting to find the closure by iterating on ℒ_β
generates a splatter of step functions. This case is examined more
closely in the next chapter. 

Attempting some guess-work, the self-similarity of the cpr function
suggests an opening. Specifically, let _β(x)=_β(x)-1/2.
The one finds that

    [𝒫_β_β](y)=    1/β[_β(y/β)+_β(y/β+1/2)]
    
    =    _β(y)/β

This is a non-polynomial, fractal eigenfunction of 𝒫_β,
and, with a bit of elbow-grease, one can find many more. This includes
the Takagi functions, and their higher-order analogs, which are, roughly
speaking, Takagi functions constructed from polynomials. These all
have interesting self-similarity properties under the dyadic monoid.

Unfortunately, one has that _β(x)0
when β<2x; it won't do as an eigenfunction of ℒ_β.
There is no obvious, simple modification of _β(x)
that would cause it to be a valid eigensolution of ℒ_β.
Manually installing a factor of Θ(β/2-y)
and then iterating to find the closure leads to the same splatter
of step functions as in the case of the polynomials.

Another interesting case arises if one attempts a Fourier-inspired
basis. Define 

    _β;n;l(x)=exp i2π(2l+1)β^nx

for integer l. One then obtains a shift sequence

    [𝒫_β_β;n;l](x)=1/β_β;n-1;l(x)(1+_β;n;l(1/2))

This is not a viable candidate for ℒ_β, as it is
again beset by the step function. As a shift sequence, it can be used
to construct coherent states that are eigenfunctions of 𝒫_β,
having any eigenvalue within the unit disk. Specifically, observe
that _β;0;l(1/2)=exp iπ(2l+1)=-1
so that [𝒫_β_β;0;l](x)=0
and so the shift sequence terminates after finite iteration. Given
a complex value z, construct the coherent state as

    ϕ_l;z(x)=∑_n=0^∞z^n_β;n;l(x)

The shift is then

    [𝒫_βϕ_l;z](x)=z/β∑_n=0^∞z^n(1+_β;n+1;l(1/2))_β;n;l(x)

This is not particularly useful, until one notices that for for certain
values of β, this contains nilpotent sub-series. 

Specifically, fix a value of n=N and consider those values of β
for which _β;N;l(1/2)=-1. This holds whenever
β^N is an odd integer, that is, whenever β=(2m+1)^1/N
(and, as always, β≤2). For these special values of β,
one has that [𝒫_β_β;N;l](x)=0
and so the functions 

    ϕ_l;z;N(x)=∑_n=0^Nz^n_β;n;l(x)

vanish after N iterations of 𝒫_β. That is, these
can be used to form a a basis in which 𝒫_β is nilpotent.
Conversely, letting m and N be free, the values for which β=(2m+1)^1/N
are dense in the interval [1,2) and so any β
is arbitrarily close to one with a nilpotent function space. These
values of β are exactly the same values for which the bit sequences
given by eqn <ref> eventually terminate in all zeros;
i.e. become periodic fixed points with period 1. 

The existence of a dense set of fixed points is dual to the the existence
of nilpotent densities. That is, one “causes” or “forces”
the other to happen. This idea should be further elaborated, as it
establishes a duality between countable and uncountable sets, which
has an element of curiosity to it.

Presumably, there are special values of β which allow a periodic
orbits to originate from a dense set. Such values of β, and
such periodic orbits, should then correspond to specific self-similarities
of the ϕ_l;z(x) function, specifically manifesting
as cyclic behavior in (1+_β;n+1;l(1/2))^p
for some period p. Whether there is some similar manifestation
for ℒ_β is wholly unclear; however, the examination
of the periodic orbits of the beta shift, undertaken in a later chapter,
will provide a strong clue.



 §.§ Rotated Renyi-Parry function


A different class of almost-solutions starts with the Renyi-Parry
invariant measure defined in eqn <ref>. Using
this, one may define a “rotated” function

    v_β;z(x)=∑_n=0^∞d_n(x)/β^nz^n

where d_n(x)=ε_n(2x/β)
are the same digits as defined by Renyi and Parry, just rescaled for
the beta shift convention being used here. That is, 

    d_n(x)=
    1    x<T^n(β/2)
    
    0

and T the beta shift map of eqn <ref>. The summation
is clearly convergent (and holomorphic) for complex numbers z within
the disk |z|<β.

Exploring this numerically, one finds that 

    [ℒ_βv_β;z](x)=1/zv_β;z(x)+C(β;z)

where C(β;z) is a constant independent of x (or
rather, a constant for x<β/2 and zero otherwise). The fact
that C(β;z) is a constant independent of x can
be taken either as a humdrum fact, typical of the territory, or as
a big unexplained surprise! I prefer the later. The sequence of digits
d_n(x) depend strongly on x. The function v_β;z(x)
depends strongly on x. The digits d_n(x) are entangled
in an analytic series, in very highly non-trivial form. This is worth
pursuing to great extent.

First, some baby-steps. Combining the above expressions, one can write

    C(β;z)=    [ℒ_βv_β;z](y)-v_β;z(y)/z
    
    =    ∑_n=0^∞z^n/β^n[1/βd_n(y/β)+1/βd_n(y/β+1/2)-1/zd_n(y)]

Again: this is a constant, independent of y for 0≤ y<β/2,
and zero for y≥β/2. This only holds for the sum; individual
terms are not independent of y. The interesting limit is where
|z|→β and so its convenient to define

    D(β;ζ)=    ζβ C(β;ζβ)
    
    =    ∑_n=0^∞ζ^n[ζ d_n(y/β)+ζ d_n(y/β+1/2)-d_n(y)]

This is holomorphic on the unit disk |ζ|<1, as each
individual d_n is either zero or one; there won't be any poles
inside the unit disk. Again, this is independent of y. 

The equation is readily simplified. Set y=0, to obtain d_n(0)=1.
Cancelling terms, one obtains

    D(β;ζ)=-1+ζ∑_n=0^∞ζ^nd_n(1/2)

This is straight-forward to examine numerically; the disk is readily
plotted, and reveals what looks like a countable number of zeros in
the disk |ζ|<1, and then many more in the limit
|ζ|→1. Each zero corresponds to an eigenfunction/eigenvalue
pair, the eigenfunction given by eqn <ref>.
The table below lists some examples; the numbers are accurate to about
the last decimal place.






β     z     |z|     1/z1.8     -1.591567859     1.591567859     -0.62831125581.8     -1.1962384 +i 1.216022231     1.705783215     -0.4111213835 - i 0.4179206604 1.8     0.9919147363 +i 1.446092984     1.753590535     0.3225655308 -i 0.47026194291.6     -1.063651387 +i 1.008959895     1.466067646     -0.4948701876 -i 0.46942464291.4     0.550836432 +i 1.178171082     1.300579822     0.3256481633 -i 0.69652119311.2     0.9578845659 +i 0.6073301155     1.134192537     0.7446284155 -i -0.4721187476





These are not particularly meaningful numbers; they just give a flavor
for some locations of eigenvalues. The overall distribution is scattered;
the zeros appear in not very predictable locations, mostly not far
from the edge of the disk. Insight into this behavior is developed
in depth in later chapters. A typical eigenfunction is shown in figure
<ref>; this is for the zero listed in the
last row of the table above. Although its unlike figure <ref>,
it does have the same general characteristics, belonging to the same
family.



Formula <ref> will reappear in a later chapter,
where it is found to describe periodic orbits of the beta shift. It
will serve to tie together the transfer operator with the behavior
of orbits.

This section ends with a deep mystery: why is C(β;z)
independent of the y coordinate? How should it be explained? How
can it be put to good use? Might there be more of such identities?
Or is there only this one? 



 §.§ Analytic Ergodics


Why is C(β;z) independent of the y coordinate?
And is this surprising, or not surprising?

Perhaps it should not be surprising. The full shift (the Bernoulli
shift) has an absolutely flat, uniform distribution. This is not a
surprise: real numbers are uniformly distributed on the real number
line. Pick almost any real number: the binary digits are uniformly
distributed. The FP eigenfunction of the transfer operator for the
Bernoulli shift is a constant. This is not a surprise; the uniformity
of the reals demands this.

Compare to the invariant measure of eqn <ref>,
shown in the very first figure <ref>.
The figure clearly consists of flat plateaus; the flatness harking
back to, and “explained by” the flatness of the real numbers.
But those flat regions arise from eqn <ref>
(with z=1), and if one looks at the individual binary digit sequences
d_n(x) for two different values of x, these sequences
become wildly, crazily, chaotically different; why should they sum
to the same value? This now seems to be profound, in some way that
the flatness of the reals is not. 

Presumably, this should be known as “the fundamental theorem of
analytic ergodics”. Or at least, “the funamental theorem of analytic
subshifts”. Oddly, I've never heard of such a theorem, nor recollect
any proof of it having scrolled past my eyes. Perhaps it is known.
It certainly deserves greater fame; it seems like it would be quite
the tool to unlock many of the mysteries being fumbled-over in this
text.



 §.§ Iterated transfer operator


To understand the nature of the steady-state solution (the Frobenius-Perron
eigenstate), its is worth iterating on the recurrence relation for
it, by hand, the first few times. To do this, it is convenient to
write it in the form

    [ℒ_βf](y)=Θ(y)/β[f(α(y))+f(ω(y))]

where Θ(y)=1 if y≤β/2 else zero; this
is a step function to denote the vanishing for the operator for 2y>β.
(This differs from the use of Θ as the Heaviside step function
in earlier sections; the intent is the same, but the goal is to have
a briefer notation here. Which is which should be clear from context.)
The functions α(y)=y/β and ω(y)=1/2+α(y)
are convenient shorthands for symbolic iteration.

Iterating once gives

    [ℒ^2f](y)=    Θ(y)/β^2   [Θ(α(y))[f(α^2(y))+f((ω∘α)(y))]+.
          .Θ(ω(y))[f((α∘ω)(y))+f(ω^2(y))]]

Using a simplified notation g(y)=f(α(y))+f(ω(y))
allows this to be iterated a third time:

    [ℒ^3f](y)=    Θ(y)/β^3   [Θ(α(y))[Θ(α^2(y))g(α^2(y))+Θ(ωα(y))g(ωα(y))]+.
          .Θ(ω(y))[Θ(αω(y))g(αω(y))+Θ(ω^2(y))g(ω^2(y))]]

and a fourth time, this time omitting the argument, and the various
nesting parenthesis.

    [ℒ^4f](y)=    Θ(y)/β^4   [ΘαΘα^2[Θα^3gα^3+Θωα^2gωα^2]+.
          ΘαΘωα[Θαωα gαωα+Θω^2α gω^2α]
          ΘωΘαω[Θα^2ω gα^2ω+Θωαω gωαω]
          .ΘωΘω^2[Θαω^2gαω^2+Θω^3gω^3]]

Notice that the primary structure is given by a product of step functions.
This is more conveniently visualized as a tree:




    < g r a p h i c s >




For any given iteration, the result is the sum of the vertexes at
a given level, while the product of step functions is the product
of the step functions in the tree, following the path to each node.
Because any particular step function might be zero, it effectively
acts to cut off the tree at that location. It is therefore interesting
to understand general products of the α and β functions. 

It is convenient to define 

    γ_x(y)=x/2+y/β

so that α(y)=γ_0(y) and ω(y)=γ_1(y),
so that a general iterated sequence of intermixed α's and
ω's can be written uniformly in terms of γ. Given
a sequence of bits b_0b_1b_2⋯ b_n with each b_k
being either zero or one, the iterated sequence of functions can be
written as

    (γ_b_0γ_b_1γ_b_2⋯γ_b_n)(y)=1/2[b_0+b_1/β+b_2/β^2+⋯+b_n/β^n]+y/β^n+1

So, for example: 

    α^n(y)=y/β^n

while 

    ω^2(y)=1/2+1/β(1/2+y/β)

and, in general, that

    ω^n(y)=1/2[1+1/β+1/β^2+⋯1/β^n-1]+y/β^n

Iterated sequences of pairs of functions, of the form γ_b_0γ_b_1γ_b_2⋯γ_b_n
are reminiscent of de Rham curves, which generalize Cesaro curves
and the Koch snowflake. The proper definition of a de Rham curve assumes
the sequence is of infinite length, and applies a certain continuity
condition, and is generally carried out on the complex plane, so that
a continuous, nowhere-differentiable curve results. Here, the curve
is distinctly not continuous: eqn <ref> is a
finite-length form of the shift series <ref> which
can be visualized as the expander function pdr <ref>,
as shown in figure <ref>. 



 §.§ The Tree Function


Given a bit sequence (b_k) and value for y, define
the tree function as 

    T_β((b_k);y)=Θ(y)∏_n=0^∞Θ(γ_b_0γ_b_1γ_b_2⋯γ_b_n(y))

For any given fixed sequence of bits and value of y, this function
is either zero or one. One way to understand this function is to ask
how it varies for fixed β and y, but with the bit sequence
coming from the Bernoulli shift of eqn <ref>, so
that b_n=b_n(x). This simplifies notation, so that
one can write

    T_β(x;y)=T_β((b_k(x));y)=Θ(y)∏_n=0^∞Θ(γ_x;n(y))

with γ_x;n(y)=γ_b_0γ_b_1γ_b_2⋯γ_b_n(y).
Its clear that the tree function has maximum support when y=0.
Figure <ref> shows several gamma functions, and
the corresponding tree function that results. Figure <ref>
shows the x vs. y behavior of the tree functions. Figure <ref>
shows the β vs. x behavior of the functions. Figure <ref>
shows a unified visualization of the three preceding charts.












 §.§ Haar Basis Matrix Elements


The symmetric Haar wavelets are built from the mother wavelet

    h(x)=
    1    0≤ x<1/2
    
    -1    1/2≤ x<1

and has individual wavelets given by 

    h_nj(x)=2^n/2h(2^nx-j) 0≤ j≤2^n-1

The matrix elements of the transfer operator are

    ⟨ mi| L_β|nj⟩ =∫_0^1h_mi(x)[ℒ_βh_nj](x)dx

where the operator ℒ_β is given by eqn <ref>.
Computing these by hand promptly pushes into a big mess. One can obtain
explicit expressions, just that they are tedious to obtain. Some preliminary
observations include that 

    ⟨ mi| L_β|nj⟩ =0 β≤ i/2^m-1

because the transfer operator vanishes above β/2. In the same
vein, matrix elements vanish unless 

    [i/2^m,i+1/2^m]∩[β j/2^n,β(j+1)/2^n]∅

or if 

    [i/2^m,i+1/2^m]∩[β(j/2^n-1/2),β(j+1/2^n-1/2)]∅

In all other cases, the Haar wavelets completely fail to overlap,
and thus the matrix elements are zero. In addition, only three pairs
of wavelets overlap in a non-zero fashion. That is, for a fixed m,n
and j, there are at most six different values of i for which
the matrix elements are non-vanishing: the first three of these are
the values for which

    β j/2^n∈[i/2^m,i+1/2^m]  β(j+1/2)/2^n∈[i/2^m,i+1/2^m]  β(j+1)/2^n∈[i/2^m,i+1/2^m]

and likewise for three more. The observation is that the integral
vanishes unless the first wavelet intersects an edge transition of
the second wavelet.

The primary failure of this basis is that there is no obvious way
to diagonalize the transfer operator in this basis. There is no obvious
way of solving it, of finding it's eigenfunctions and eigenvalues,
other than by brute-force numerical attack.



 §.§ Julia Set


Consider the two iterators a_0(y)=min(β/2,β y)
and a_1(y)=max(0,β y-β/2).
Individually, they are the two arms of the beta shift. Here, they
have been separated from each other, so that the full domain 0≤ y≤1
is allowed. Exploring all possible inter-iterations for these gives
the Julia set for the transfer operator: it indicates where a point
“came from”, for the iterated transfer operator. There are several
related ways to visualize this. One way is to fix y and then, given
a bit-sequence (b_n) to compute 

    j((b_n))=a_b_o∘ a_b_1∘ a_b_2∘⋯(y)

Figure <ref> shows a visualization for
finite bit-sequences: in essence, the very first few iterations. Although
it is similar to figure <ref>, it is not the same. 



For a related notion, consider the definition of “laps”, from
Jeffrey Lagerias etal.<cit.>. 



§ HESSENBERG BASIS


There is a set of Haar-like wavelets in which the transfer operator
is of the form of a Hessenberg operator - that is, the operator becomes
almost upper-diagonal, with only one diagonal, just below the main
diagonal, that is non-zero. Explicitly, the transfer operator ℒ_β
has matrix entries [ℒ_β]_ij such that
[ℒ_β]_ij=0 whenever i>j+1. A matrix
having this form is called a Hessenberg matrix; such matrices have
various interesting properties; among others, they generalize the
Jacobi matrix. This chapter explicitly constructs an infinite-dimensional
Hessenberg matrix, which may now be called a Hessenberg operator.

Hessenberg operators occur naturally in spectral measure theory; some
of this will be reviewed in several later chapters. To get a flavor
for what is to come: Given a Hessenberg operator, one may construct
a system of orthogonal polynomials that provide a basis for square-integrable
holomorphic functions on some domain of the complex plane. Such a
space is called a Berman space; in this sense it generalizes the Jacobi
operator for real Borel measures. This basis of polynomials in turn
allows the Hessenberg operator to be explicitly seen as a shift operator
on that domain, with [Hf](z)=zf(z)
for H the Hessenberg operator and f(z) a holomorphic
function (specifically a Bergman function) on the Bergman domain.
But all of this is for later chapters; its mentioned here only to
whet the appetite. 



 §.§ Hessenberg wavelet basis


The transfer operator ℒ_β can be fairly easily
brought into Hessenberg matrix form. A sequence of of orthonormal
functions is constructed in this section; when used as a basis, the
transfer operator becomes almost upper-diagonal. 

The trick to the construction is to define wavelets such that the
transfer operator applied each wavelet causes the end-points of the
wavelet to exactly line up with the end- or mid-points of previous
wavelets, thus avoiding the nasty interval-overlap algebra required
with the Haar basis. This is accomplished by carefully picking the
midpoint of the next wavelet in the sequence to be located exactly
at the discontinuity of the transfer operator applied to the previous
wavelet.

The construction proceeds as follows. Let 

    ψ_0(x)=1/√(β/2)   0≤ x≤β/2
    
    0    β/2<x≤1

Consider L_βψ_0. It is the sum of two parts:
two step-functions; one which is constant for x≤β/2 and another
that is constant for x/β+1/2≤β/2.
Solving explicitly for the location of the step, it is x=β(β-1)/2.
For convenience, define m_1=β(β-1)/2 and m_0=β/2.
These will anchor a series of midpoints, beginning with m_-1=0.
Using the midpoint m_1, construct the wavelet

    ψ_1(x)=1/m_1√(m_1(m_0-m_1)/m_0)   0≤ x≤ m_1
    -1/m_0-m_1√(m_1(m_0-m_1)/m_0)   m_1<x≤ m_0
    
    0    m_0<x≤1

Note that this is normalized to unit length: ∫_0^1|ψ_1(x)|^2dx=1
and that it is explicitly orthogonal to the first: ∫_0^1ψ_1(x)ψ_0(x)dx=0.

Consider L_βψ_1. As always, it is the sum of
two parts. The midpoint of ψ_1 is at m_1=β(β-1)/2
and this mid-point is mapped to one of two different places. If m_1<1/2
then it is mapped to m_2=β m_1 else it maps to m_2=β(m_1-1/2).
Thus, if m_1<1/2, define 

    ψ_2(x)=
    0    0≤ x≤ m_1
    1/(m_2-m_1)√((m_2-m_1)(m_0-m_2)/m_0-m_1)   m_1≤ x≤ m_2
    -1/(m_0-m_2)√((m_2-m_1)(m_0-m_2)/m_0-m_1)   m_2<x≤ m_0
    
    0    m_0<x≤1

else define 

    ψ_2(x)=1/m_2√(m_2(m_2-m_1)/m_1)   0≤ x≤ m_2
    -1/(m_1-m_2)√(m_2(m_2-m_1)/m_1)   m_2≤ x≤ m_1
    
    0    m_1<x≤1

Because each end of the interval on which ψ_2 is non-zero
lies entirely within one of the constant arms of ψ_1, one
has, by construction, that ∫_0^1ψ_2(x)ψ_1(x)dx=0
(and, of course, ∫_0^1ψ_2(x)ψ_0(x)dx=0.)

The rest of the basis can be constructed iteratively, based on these
examples. The midpoints are given by iterating <ref>
on m_0=β/2, so that m_p=T_β(m_p-1)=T_β^p(m_0)
is the p'th iterate of β/2. Let m_l be largest midpoint
smaller than m_p (and l<p); let m_u be the smallest midpoint
larger than m_p (and l<p). Let m_-1=0 initiate the sequence
by providing the smallest-possible “midpoint”; m_0=β/2
already provides the largest possible. 

Then define 

    ψ_p(x)=
    0    0≤ x≤ m_l
    C_p/(m_p-m_l)   m_l≤ x≤ m_p
    -C_p/(m_u-m_p)   m_p<x≤ m_u
    
    0    m_u<x≤1

By construction, this has the property that ∫_0^1ψ_p+1(x)ψ_n(x)dx=0
for any n<p+1. The normalization constant is 

    C_p=√((m_p-m_l)(m_u-m_p)/m_u-m_l)

which is determined by requiring that ∫_0^1|ψ_p(x)|^2dx=1.



 §.§ Matrix Elements


The above-defined basis provides the Hessenberg representation for
the transfer operator. Defining

    ⟨ n| L_β|m⟩ =∫_0^1ψ_n(x)[ℒ_βψ_m](x)dx

this has the expected Hessenberg form, in that 

    ⟨ n| L_β|m⟩ =0  n>m+1

This is just one diagonal short of being actually solvable. A visualization
of the matrix elements is shown in figure <ref>.






 §.§ Completeness


The Hessenberg basis construction gives a countable set of ψ_n
that is an orthonormal basis on the unit interval: ∫_0^1ψ_m(x)ψ_n(x)dx=δ_mn.
Are they complete? Obviously the {ψ_n} cannot
be complete on the unit interval, as they all vanish for β/2<x.
Perhaps they are complete on the interval [0,β/2], where they
are already orthonormal: ∫_0^β/2ψ_m(x)ψ_n(x)dx=δ_mn.

A numerical exploration shows that the midpoints m_p are dense
in the interval (0,β/2), and so this suggests that the basis
should be considered to be “sufficiently complete” on the interval
[0,β/2]. The distribution of the m_p follow exactly the
distribution of the invariant measure. Convergence is uniform to the
same degree that the midpoints “eventually” fill in and become
dense in some interval. Renyi<cit.> and Parry<cit.>
do more: they show that the midpoint process is ergodic (Parry points
out that it's weakly mixing), and provide a formal proof that the
distribution is one and the same as the invariant measure.

The above has some exceptions: there are some values of β for
which the midpoint m_0 iterates x=1/2, wherupon iteration
stops (i.e. iterates to zero), or becomes cyclic (forming a periodic
orbit). Which is which depends on how the point 1/2 is treated
by the map. These values of β are potential “trouble spots”,
and are explored in greater detail in the next chapter. They are dense
in the interval 1<β<2, but they form a countable set that can
be taken to be of measure zero. Thus, most “most” values of β
are not problematic. Excluding the trouble spots, the Hessenberg basis
can be taken to be complete.

Clearly, the ψ_n span some subspace; do they span the Hilbert
space L_2[0,β/2] of square-integrable functions
on the interval [0,β/2]? To what degree can one legitimately
write

    δ(y-x)=∑_n=0^∞ψ_n(y)ψ_n(x)

as a resolution of the identity? 

The question of completeness dogs some “obvious” assumptions one
wants to make. For example, if the set of states is complete, and
the resolution of the identity holds, then one expects that the transfer
operator resolves to the iterated function:

    δ(y-(β x 1))=∑_n=0^∞∑_m=0^∞ψ_n(y)⟨ n| L_β|m⟩ψ_m(x)

It is fun to verify that the world works as one expects it to work:
the above can be verified to hold numerically, for sums limited to
a finite cutoff.



 §.§ Numerical Eigenvalues


Given the apparent sparsity visible in figure <ref>,
one might think that the eigenvalue problem is fairly stable, numerically.
It is not all that much. Numerical exploration suggests that the spectrum
is on or near a circle lying in the complex plane[This was confirmed with both GSL gsl_eigen_nonsymmv() and Lapack
DHSEQR solvers, exploring the principle submatrix of various sizes,
up to about 2000×2000 entries. Both systems agree to at least
six decimal places, if not more. Both show sporadic eigenvalues off
the circle, but these are not numerically stable; ergo, the only valid
eigenvalues are those on the circle. The matrix entries were constructed
using the midpoint algorithm, described in the last section. To verify
that they are accurate, several techniques were used: numerical integration
to confirm orthogonality, and the use of the GnuMP multi-precision
library to push up accuracy.], of radius |λ|=1/β (ignoring, that is, the
leading eigenvalue of 1, which is easily found). 

To be clear, this is a numerical exploration of the N× N principle
submatrix of ⟨ n| L_β|m⟩.
The eigenvalue problem being posed is to find a vector v⃗=(v_k)_k=0^N
that solves 

    ∑_m=0^N⟨ n| L_β|m⟩ v_m=λ v_n

for some constant λ (with the set of possible λ
depending on N, of course). 

There are various pitfalls in extrapolating from this to the N→∞
limit. For the next few paragraphs, consider only some notion of a
“minimal” extension from finite N to the limit. That is, for
each finite N, one has a finite set of eigenvalues and eigenvectors.
In the limit, there may be accumulation points: points where the eigenvalues
accumulate to a limit point, in a standard topological sense. What
should that that topological space be? For finite N, all eigenvectors
are explicitly summable, and thus can be taken to belong to any Banach
space ℓ_p. One may as well take p=2 the Hilbert space,
and normalize the eigenvectors v⃗ so that 1=∑_m=0^Nv_m^2. 

For finite N, it appears that “most” eigenvalues λ
are “near” the circle |λ|=1/β, and that
they seem to be very uniformly distributed around this circle. The
numerical results indicate that in the limit N→∞, that the
scare-quotes “most” becomes “almost all” in the usual sense.
Similarly, “near” appears to mean that for any given λ
at finite N, one has that |λ|-1/β∼𝒪(1/N).
As to uniformity, it seems that the spacing between nearest neighbors
is also 𝒪(1/N), and that there are no “premature”
accumulation points: eigenvalues never get any closer than 𝒪(1/N),
either.

Thus, the minimal closure, the minimal extrapolation to limit points
strongly suggests that the limit points really do lie, uniformly distributed,
on the circle |λ|=1/β. Then, writing a given
accumulation point as λ=β^-1exp2π iϕ, what the
numerics do not reveal, or, at least, do not easily reveal, is whether
the allowed values of ϕ are always rational, irrational or might
have some other regular structure. The numerical exploration does
suggest that the eigenvalues are dense on the circle. Certainly it
is the case Hessenberg basis is countable, an so one would expect
the eigenvalue spectrum obtained in this way to be at least countable,
as well. Whether it is also uncountable seems unknowable in this naive
sense.

This question is interesting because if only rational ϕ are
allowed, then the decaying eigenfunctions belong to a cyclic group,
and exhibit an exact form of Poincaré recurrence as they decay. If
irrational ϕ are allowed, then the decaying eigenfunctions are
at least ergodic.

For β=2, the β-transform is the Bernoulli shift, the
transfer operator is solvable, and the spectrum is exactly known.
This has been explored by various authors<cit.>. I've written
extensively about this spectrum and the eigenvalues in other texts<cit.>.
To recap, it takes several forms, depending on the function space
that one chooses to work in. If one restricts oneself to polynomial
eigenfunctions, then the spectrum is real, non-negative (it has an
extensive kernel) and has eigenvalues of 2^-n for all n. The
eigenfunctions are the Bernoulli polynomials. Restricting to square-integrable
eigenfunctions, the spectrum continuous, having eigenvalues on the
unit disk in the complex plane. The continuous-spectrum eigenfunctions
(for eigenvalues other than 2^-n) can be understood in several
ways: if forced to be differentiable, then they are not bounded (they
diverge) at the endpoints of the interval. If forced to be bounded,
then they are fractal (non-smooth) over the entire interval. The unitary
spectrum corresponds to differentiable-nowhere eigenfunctions (wait,
or continuous-nowhere? I forget.)

A pair of plausible, numerically-extracted eigenfunctions are shown
in image <ref>.



Presumably, the spectrum can be related to the lap-counting function,
given by Lagarias<cit.>.



 §.§ (Non-)Unitarity


The numerical results suggest a hypothesis that perhaps some fragment
of ℒ_β is unitary, as it is ordinarily the case
that when eigenvalues appear on the unit circle, it is because an
operator is unitary. That does not seem to be the case here. Specifically,
define the Frobenius-Perron eigenvector ρ as the one satisfying
ℒ_βρ=ρ and normalizing it to unit length,
so that ‖ρ‖ =1 in the Hilbert (mean-square)
norm. Define the reduced operator ℛ_β in terms
of the matrix elements

    1/β⟨ n| R_β|m⟩ =⟨ n| L_β|m⟩ -⟨ρ| n⟩⟨ρ| m⟩

That is, it is just the beta shift operator, with the Frobenius-Perron
eigenvector removed, so that ℛ_βρ=0 . Its rescaled,
so that the remaining eigenvectors of ℛ_β lie on
the unit circle. Is this operator unitary in any way? That is, might
either ℛ_βℛ_β^† or ℛ_β^†ℛ_β
be the identity? Here, the dagger † is just the transpose,
as ℛ_β is purely real. Numerical exploration clearly
shows that ℛ_β is neither unitary on the left nor
on the right. Not a surprise, but does leave the understanding of
ℒ_β in a curious situation.

Perhaps it is not enough to subtract the invariant measure: The zeros
of the formula <ref> lying inside the disk must
be subtracted as well. There seems to be a cutable number of these;
the subtraction won't be straight-forward.



 §.§ Invariant Measure


Let v_n be the Ruelle-Frobenius-Perron eigenvector in the Hessenberg
basis. That is, let v_n be the vector that solves

    ∑_m=0^∞⟨ n| L_β|m⟩ v_m=v_n

This is readily computed numerically, and it is straightforward to
verify the numerics by confirming that 

    ρ(x)=∑_m=0^∞v_mψ_m(x)

is the invariant measure of equations <ref>,<ref>,
with the ψ_k(x) being the wavelets of eqn <ref>.
This expansion seems to “make sense”, as the discontinuities seen
in the graph of ρ(x) in figure <ref>
occur at exactly the midpoints m_p and the size of each discontinuity
appears to get smaller as p gets larger. Given that the wavelet
ψ_p(x) has its central discontinuity at m_p
and is bounded on left and right by midpoints of lower order, this
expansion seems to be very natural. This is supported by the diagram
<ref>, which depicts the values of v_n as a
function of n for selected values of β. These values of
v_n are real, positive, and quickly get small; there are no difficulties
or issues of convergence.



Is there some simple expression for the values of v_n as a function
of β? If so, it must be formed using some sort of fractal shift.
Figure <ref> illustrates v_1 through v_5.



The orbit of the midpoint is correlated with value of the coefficients,
illustrated in figure <ref>. The midpoint
polynomial for m_p=T_β^p(β/2), given in
eqn <ref>, is compared to v_mβ^m. It
can be seen to “line up”. The two are somehow related; its not
clear just how.






 §.§ Generating Function


The truncated ordinary generating function associated with the eigenvector
of eqn <ref> is

    G_N(z)=∑_m=0^Nv_mz^m

with the ordinary generating function being the limit N→∞.
A numerical study of this function indicates that most of the N
zeros of G_N are arranged approximately on a circle of radius
β. The arrangement appear to be quite uniform, with more-or-less
equi-distant spacing of the zeros. As N increases, it seems that
more of the zeros get closer to the circle, although the numerical
instabilities associated with double-precision math make this difficult
to control; an arbitrary-precision eigenvalue solver would be needed
to confirm this behavior.

If this behavior persists, and it seems that it will, then the limit
N→∞ cannot be taken, and the ordinary generating function,
as an analytic function, can't exist, per se, as it would be
uniformly zero inside the disk. Thus, the zeros already found by means
of eqn <ref> seem to come to the rescue: these
are located inside the disk; perhaps these are masquerading as “numerical
instabilities”, and should be taken as actually existing, and not
spurious.

In the next chapter, it will be seen that circles of zeros in the
complex plane is a recurring theme. This suggests a hypothesis that
somehow it might hold that

    ∑_m^Nv_m(β z)^m∼ z^N+1-∑_kb_kz^k

as both sides have zeros arranged in circles of unit radius. The right
hand side is defined and explored in detail in the next chapter. Superficially,
this hypothesis is clearly false: coefficients on the left are all
real and positive; coefficients on the right - the b_k, are bits,
either zero or one. Yet both exhibit a circle of zeros.

XXX This section is awkward. Revise it or cut it.



 §.§ Givens rotations


An open question: A Hessenberg matrix can be brought to solvable form
by applying a sequence of Givens rotations. Is the sequence of angles
that appear in these rotations meaningful in any way, or are they
just some form of uninteresting junk?

 



§ PERIODIC ORBITS


The iteration of the midpoint, used to construct the Hessenberg basis,
works well, unless the midpoint iterates to hit the point x=1/2
where the map has a discontinuity. Here, iteration stops: at the next
step, this point is defined to iterate to zero, in eqn <ref>.
Zero is a fixed point, and so there is nowhere further to go. This
occurs for certain values of β: after a finite number for steps,
the midpoint m_0=β/2 iterates to 1/2. This section explores
these special values of β.

Aside from the definition in eqn <ref>, one can consider
the modified map, where the less-than sign has been altered to a less-than-or-equals:

    T_β^≤(x)=β x    0≤ x≤1/2
    β(x-1/2)    1/2<x≤1

In this map, the point x=1/2 iterates to β/2, which is just
the initial midpoint itself. In this case, the halted orbits become
periodic orbits. There is a third possibility, to simply remove the
points 0, 1 and 1/2 from the domain:

    T_β^<(x)=β x    0<x<1/2
    β(x-1/2)    1/2<x<1

In this case, if the midpoint iterates to 1/2, it can be taken to
simply have exited the domain of validity.

All three variants can be considered together, so that the “true”
beta shift is taken as the quotient space or identification space<cit.>
of the three variants, in the strict topological sense of a quotient
space. Thus, interestingly, for the beta shift, the periodic orbits
and the fixed point both belong to the same equivalence class. This
has some interesting implications when one compares the beta shift
to other iterated maps, such as the logistic map, which have non-trivial
stable regions. Topologically, it would seem that one can perform
a kind of surgery, attaching stable regions exactly into those spots
where, in the beta shift, one has an equivalence class. This solves
(at least for me) the long-standing problem of exactly how to properly
describe the topological conjugacy between different kinds of iterated
maps.



 §.§ The β-generalized Golden Ratio


The above wavelet basis seems to be well-behaved, except when β=φ=(1+√(5))/2
the Golden Ratio. In this situation, one has that m_0=φ/2
and m_1=1/2. At this location, further iteration breaks down.
That is, m_2=T_φ(m_1) can either be taken
to be m_2=0 or m_2=m_0. In the former case, iteration
simply stops; in the later case, it repeats, again, without generating
new midpoints that can provide a workable basis.

Working backwards, this issue re-appears whenever the p'th iterate
m_p=T_β^p(m_0) lands at the discontinuity,
so that one may take either m_p=0 or m_p=m_0. For p=3,
there are two trouble spots, which occur when either β^3-β^2-1=0
or when β^3-β^2-β-1=0. These correspond to the
values of β=1.465571231876768⋯ and β=1.839286755214161⋯.

Where are the trouble spots located? Consider, for example, m_4=T_β^4(m_0),
and consider the movement of m_4 as β is swept through
the range 1<β<2. This is shown in figure <ref>.
As made clear in the image, three new trouble spots appear. These
are located at β=1.380327757⋯ and β=1.754877666⋯
and β=1.927561975⋯, which are the real roots of β^4-β^3-1=0
and β^4-β^3-β^2-1=0 and β^4-β^3-β^2-β-1=0
respectively.



Following a similar suggestion by Dajani<cit.>, numbers
of this kind may be called “generalized golden means”. Unfortunately,
the term “generalized golden mean” is in common use, and is applied
to a variety of different systems. Not all are relevant; one that
is, is given by Hare et al.<cit.> who provide series
expansions for the real roots of β^p-∑_k=0^n-1β^k=0;
these are known as the n-bonacci constants (Fibonacci, tribonacci,
tetranacci, etc.). Stakhov<cit.> considers β^p+1-β^p-1=0
in general settings. Some, but not all of these numbers are known
to be Pisot numbers or Salem numbers<cit.>. In what follows,
these will be referred to as the “beta golden means”, since all
of the ones that appear here have explicit origins with the beta shift.



 §.§ Counting Orbits


How many trouble spots are there? The table below shows the count
M_p of the number of “new” trouble spots, as a function of
the midpoint index p. 



p     2     3     4     5     6     7     8     9     10     11     12M_p     1     2     3     6     9     18     30     56     99     186     335



This appears to be Sloane's OEIS A001037 which has a number of known
relationships to roots of unity, Lyndon words, and the number of orbits
in the tent map. The trouble spots are the positive real roots of
polynomials of the form β^p-β^p-1-b_1β^p-2-b_2β^p-3-⋯-1=0
with the b_k being certain binary bit sequences. There is just
one such (positive , real) root for each such polynomial. These polynomials
are irreducible, in the sense that a bit-sequence b_k is disallowed
if it has the same root as some lower-order polynomial. For example,
β^4-β^3-β-1 is disallowed; it has the same root
as β^2-β-1. Although the digits b_k must be zero
or one, this definition of irreducibility, plus the counting, suggests
some relationship to the irreducible polynomials over the field _2,
as that is what the definition of OEIS A001037 counts. Yet the relationship,
if any, is quite unclear.[A hypothesis presented in a later section suggests that each orbit
should be thought of as a Galois group, with the length of the orbit
corresponding to the number of elements in the Galois group. It seems
that this might explain much of the structure.]

The values of M_n are given explicitly by Moreau's necklace-counting
function

    M_n=1/n∑_d|n2^dμ(n/d)

where the sum runs over all integers d that divide n and μ
is the Möbius function. The generating function is

    t/1/2-t=∑_n=1^∞nM_nt^n/1-t^n

which has a radius of convergence of |t|<1/2. For large
n, the asymptotic behavior can be trivially deduced from the defining
sum:

    M_n=2^n/n-𝒪(2^n/2/n)


The above counting function is for necklaces with only two colors.
In general, one can have necklaces with 3 or more colors; can that
happen here? Yes, of course: if one considers the general β-transform
for 2<β, then, in general, it can be taken as a “kneading
transform” with ⌈β⌉ branches or folds
in it. The analogous trouble-spots again appear, and they can appear
after an arbitrary finite-length orbit. Insofar as they correspond
to periodic orbits, they are necessarily counted by the necklace-counting
function. That is, one must consider all possible strings of ⌈β⌉
letters, modulo a cyclic permutation: this is the very definition
of a necklace (or “circular word”). The number of such necklaces
is given by the necklace-counting function. Each such orbit is necessarily
represented by a Lyndon word, which is a representative of the conjugacy
class of the orbit.



 §.§ β-Golden Polynomials


The “trouble spots” whenever the k'th iterate m_k=T_β^k(m_0)
lands on the starting midpoint m_k=m_0. Because of the piece-wise
linear form of T_β, the k'th iterate will be a piece-wise
collection of polynomials, each of order k, and of the form β^k-β^k-1-b_1β^k-2-b_2β^k-3-⋯-1
for some binary digits b_i being zero or one. These must be arranged
in the manner such that β^k-β^k-1-b_1β^k-2-b_2β^k-3-⋯-1=0
at each discontinuity, as illustrated in figure <ref>.
This limits the polynomials that can appear; it limits the possible
coefficients b_i; not all bit-sequences appear. 

Although the count of these polynomials is the same as that for the
irreducible polynomials over 𝔽_2, the relationship is
completely opaque (to me). Apparently, this is a generic issue: the
number of irreducible polynomials over 𝔽_p, for p
prime, is the same as the number of necklaces, and yet, there is no
known bijection between these irreducible polynomials and the Lyndon
words!<cit.> Thus, a compilation seems to be called for.

The table below explicitly shows the polynomials for the first few
orders. A polynomial is included in the table if it is an iterate
of a previous polynomial, and if it's real root is bracketed by the
roots of the earlier iterates. That is, must
have the form 
when n is even or 
for n odd. The roots must be bracketed by the roots of polynomials
occurring earlier in the sequence; if the root is not bracketed, then
the corresponding polynomial does not appear in the list.

The bracketing relationship is rather awkwardly expressed in the following
pseudo-code. Here, r_n is the root p_n(r_n)=0.
The polynomial p_n is included in the list if it is the case
that this pseudo-code does not fail:


[t]0.8















The above is a rather awkward way of stating that roots must be bracketed
by pairs of previous roots. It can perhaps be more easily understood
by studying the location of the discontinuities in figure <ref>:
new discontinuities at higher orders must occur before earlier ones.

Thus, for example, the polynomial β^3-β-1 is excluded
from the list simply because it is not an iterate of an earlier polynomial,
even though it has the interesting real root 1.324717957244746⋯,
the “silver constant”. The numbering scheme does not even have
a way of numbering this particular polynomial. Despite this, the silver
constant does appear, but a bit later, as the root of p_8=β^5-β^4-1,
which is an allowed polynomial.

The polynomial p_5=β^4-β^3-β-1 is excluded because
it has φ=1.618⋯ as a root, which was previously observed
by p_1. The polynomial p_9=β^5-β^4-β-1 is
excluded because it's root, r_9=1.497094048762796⋯ is greater
than its predecessor r_2; the recursive algorithm does not compare
to r_4. 

There are other ways to start the table, and to index the polynomials.
The given indexing has the property that 2n+1, taken as binary,
encodes the coefficients of the polynomial. The order of the polynomial
is ⌈log_2(2n+1)⌉. The index
n itself encodes the orbit of the midpoint. That is, writing n=b_0b_1b_2⋯ b_p
for binary digits b_k, then T_β^k(β/2)<1/2
if and only if b_k=0. Note that b_0=1 always corresponds
to 1/2<β/2 always. By convention, the last digit is always
1, also.




order     p_n(β)     n     binary     root0     1               2*1     β          0     02-5 
     β-1     0     1     12     β^2-β-1     1     11     φ=1+√(5)/2=1.618⋯2*3     β^3-β^2-1     2     101     1.465571231876768⋯2-5 
     β^3-β^2-β-1     3     111     1.839286755214161⋯3*4     β^4-β^3-1     4     1001     1.380277569097613⋯2-5 
     β^4-β^3-β^2-1     6     1101     1.7548776662466924⋯2-5 
     β^4-β^3-β^2-β-1     7     1111     1.9275619754829252⋯6*5     β^5-β^4-1     8     10001     1.324717957244746⋯2-5 
     β^5-β^4-β^2-1     10     10101     1.5701473121960547⋯2-5 
     β^5-β^4-β^3-1     12     11001     1.704902776041646⋯2-5 
     β^5-β^4-β^3-β-1     13     11011     1.812403619268042⋯2-5 
     β^5-β^4-β^3-β^2-1     14     11101     1.888518845484414⋯2-5 
     β^5-β^4-β^3-β^2-β-1     15     11111     1.965948236645485⋯



The next table lists the acceptable polynomials for order 5, 6 and
7. Again, the coefficients appearing in the polynomial are encoded
by the binary value of 2n+1 in the sequence.



order     sequence5     8,10,12,13,14,156     16,20,24,25,26,28,29,30,317     32,36,40,42,48,49,50,52,53,54,56,57,58,59,60,61,62,63



Although there are as many of these polynomials as there are Lyndon
words, there is no obvious way to write a bijection between the two.
It is almost possible to do so by writing 2n in binary, and then
reversing the order of the bits, left-to-right. One almost gets the
Lyndon words in the correct order, except “in the middle”: so,
for example, in the table above, one can get the Lyndon order by exchanging
10 with 12, and 13 with 14. But the table above cannot be re-ordered:
the given ordering encodes the orbit of the midpoint. Apparently,
although a given orbit can be cyclically rotated to obtain a Lyndon
word, the initial segment of the orbit is not a Lyndon word itself.

Questions that present themselves include: 


  * Is there a generating function for the squence of allowed values of
n? What is it?

  * How long is the initial segment of each periodic orbit, before the
orbit attains it's Lyndon word form? What are the values of β
where the initial orbits are not in Lyndon form?




 §.§ Distribution of β-Golden Roots


It seems natural to assume that the real roots have some distribution.
This seems not to be the case. Figure <ref>
shows the numerically computed (bin-counted) distribution of the zeros
of p_n(β) for n<2^k for three different values
of k. This suggests that, in the limit of k→∞, almost
all p_n(β) have roots that approach 2. Although
the roots appear to be dense in 1<β<2, essentially all of the
weight of that density is at 2. Since the roots are countable,
the density clearly becomes very thin.



The local distribution of roots can be sensed from the figure <ref>,
which visualizes the distance between neighboring roots.






 §.§ Complex Roots


What are the complex roots? Numerical work clearly indicates that
they seem to be approximately cyclotomic in some sense or another.
They seem to be more-or-less uniformly distributed in a circle, always.
The modulus of most of the complex roots appear to be less than one.
This is violated for the complex roots of p_2^k(β)=β^k+2-β^k+1-1,
where some of the roots in the right-hand quadrant have a modulus
larger than one. By contrast, the complex roots of p_2^k-1(β)=β^k+1-∑_j=0^kβ^j
seem to always have a modulus less than one. These two seem to be
the extreme cases: in general, the polynomials appear to be “approximately
cyclotomic”. Its not clear how to make this statement more precise.

These numerical results can be argued heuristically: just divide the
polynomial by it's leading order. That is, a general polynomial of
this form is

    p_n(z)=z^k+1-∑_j=0^kb_jz^k-j

with the convention that b_0=b_k=1, and the bit-sequence n=b_0b_1b_2⋯ b_p
corresponding to a terminating orbit. Dividing by z^k+1gives
a series 

    1-1/z-b_1/z^2-b_2/z^3-⋯

Clearly, this can have a zero only when |z|<2 as otherwise,
the terms get too small too quickly. 



 §.§ β-Golden β-Fibonacci Sequences


It is well known that the golden ratio occurs as limit of the ratio
of adjacent Fibonacci numbers:

    φ=lim_m→∞F_m/F_m-1

where F_m=F_m-1+F_m-2. There is a generalization of this,
which also has received attention: the tribonacci, quadronacci, etc.
sequences, whose limits are 

    α_n=lim_m→∞F_m^(n)/F_m-1^(n)

where F_m^(n)=F_m-1^(n)+F_m-2^(n)+⋯+F_m-n^(n). 

Is it possible that the real roots of the polynomials p_n(β)
are also the roots of such sequences? But of course they are! Given
a finite string of binary digits { b} =b_0,b_1,⋯,b_k,
write the beta-Fibonacci sequence as

    F_m^{b}=b_0F_m-1^{b}+b_1F_m-2^{b}+⋯ b_kF_m-k^{b}

The name “beta-Fibonacci” is needed because the term “generalized
Fibonacci sequence” is already in wide circulation for the special
case of all bits being one. The ratio of successive terms is 

    α^{b}=lim_m→∞F_m^{b}/F_m-1^{b}

and is given as the (positive) real root of the polynomial

    p_n(β)=β^k+1-b_0β^k-b_1β^k-1-⋯-b_k=0

These polynomials and their roots were already enumerated and tabulated
in the previous section. 

The beta-Fibonacci sequences do not appear by accident: these sequences
have an ordinary generating function (OGF) given by the polynomial!
That is,

    ∑_m=0^∞z^mF_m^{b}=z^k/1-b_0z-b_1z^2-⋯-b_kz^k+1=1/zp_n(1/z)

The factor of z^k in the numerator serves only to initiate the
sequence so that F_0^{b}=⋯=F_k-1^{b}=0 and F_k^{b}=1.

These sequences are generic: they indicate how many different ways
one can partition the integer m into elements of the set { b_0,2b_1,3b_2,⋯,(k+1)b_k}.
So, for example, the entry for n=12 in the table below corresponds
to OEIS A079971, which describes the number of ways an integer m
can be partitioned into 1, 2 and 5 (or that 5m can be partitioned
into nickels, dimes and quarters). This corresponds to the bit sequence
{ b} =11001; that is, { b_0,2b_1,3b_2,⋯,(k+1)b_k} ={ 1·1,2·1,3·0,4·0,5·1} ={ 1,2,5}.
From such partitions, it appears that one can build partitions of
the positive integers that are disjoint, and whose union is the positive
integers. This suggests a question: can these partititions be expressed
as Beatty sequences?

The previous table is partly repeated below, this time annotated with
the OEIS sequence references.






n     binary     root     root identity     sequence0     1     1          1     11     φ=1+√(5)/2=1.618⋯     golden ratio     Fibonacci A0000452     101     1.465571231876768⋯      OEIS A092526     Narayana A0009303     111     1.839286755214161⋯     tribonacci A058265      tribonacci A0000734     1001     1.380277569097613⋯     2nd pisot A086106     A0032696     1101     1.754877666246692⋯     OEIS A109134     A0609457     1111     1.927561975482925⋯     tetranacci A086088     tetranacci A0000788     10001     1.324717957244746⋯     silver A060006     A00352010     10101     1.570147312196054⋯     pisot A293506     A06096112     11001     1.704902776041646⋯          A07997113     11011     1.812403619268042⋯          A07997614     11101     1.888518845484414⋯          A07997515     11111     1.965948236645485⋯     pentanacci A103814     A001591





All of these integer sequences and roots participate in a number of
curious relations having a regular form; this is, of course, the whole
point of listing them in the OEIS. This suggests a question: do the
known relationships generalize to the beta-shift setting?

For example, there are various known relations for the “generalized
golden means”. These are the roots of the series for which all b_k=1,
that is, the roots of 

    β^k+1-β^k-β^k-1-⋯-1=0

In the present notation, these would be the roots of the polynomials
p_n(β)=0 for n=2^k-1. Such roots can be rapidly
computed by a series provided by Hare, Prodinger and Shallit<cit.>:

    1/α_k=1/2+1/2∑_j=1^∞1/jj(k+1)  j-11/2^j(k+1)

This series is obtained by making good use of the Lagrange inversion
formula. Here, α_k is the k'th generalized golden mean,
i.e. the solution p_2^k-1(α_k)=0. Can the
Hare series be extended to provide the roots r_n of p_n(r_n)=0
for general n?

Another set of observations seem to invoke the theory of complex multiplication
on elliptic curves, and pose additional questions. So:

The tribonacci root r_3 is given by

    r_3=1/3(1+√(19+3√(33))+√(19-3√(33)))≃1.839⋯


The silver number (plastic number) r_8 is given by 

    r_8=1/6(√(108+12√(69))+√(108-12√(69)))≃1.324⋯


The Narayana's cows number r_2 is given by 

    r_2=1/6√(116+12√(93))+2/3√(116+12√(93))+1/3≃1.645⋯


The root r_6 is related to the silver number r_8 as r_8=r_6(r_6-1)
and is given by

    r_6=1/6√(108+12√(69))+2/(√(108+12√(69)))^2≃1.754⋯


Do the other roots have comparable expressions? To obtain them, is
it sufficient to articulate the theory of “complex multiplication”
on elliptic curves? The appearance of only the cube and square roots
is certainly suggestive of an underlying process of points on elliptic
curves.



 §.§ β-Fibonacci sequences as shifts


The nature of the β-Fibonacci sequences as shift sequences
can be emphasized by noting that they arise from the iteration of
a (k+1)×(k+1) matrix in lower-Hessenberg
form:

    B=[[   b_0     1     0     0     ⋯     0;   b_1     0     1     0     ⋯     0;   b_2     0     0     1     ⋯     0;     ⋮     ⋮     ⋮           ⋱     ⋮; b_k-1     0     0     0     ⋯     1;   b_k     0     0     0     ⋯     0 ]]

That is, the m'th element of the sequence is obtained from the
m'th iterate B^m. That such iteration results in integer sequences
has long been observed in the theory of continued fractions. It's
useful to work an explicit example. For the golden ratio, one has

    B=[[ 1 1; 1 0 ]]

and the iterates are 

    B^2=[[ 2 1; 1 1 ]], B^3=[[ 3 2; 2 1 ]], B^4=[[ 5 3; 3 2 ]], B^n=[[   F_n F_n-1; F_n-1 F_n-2 ]]

with F_n being the n'th Fibonacci number, as usual. For the
general case, one has that

    B^m=[[     F_m^{ b}   F_m-1^{ b}   F_m-2^{ b}            ⋯ F_m-k+1^{ b}   F_m-k^{ b};   F_m-1^{ b}   F_m-2^{ b}   F_m-3^{ b}            ⋯   F_m-k^{ b} F_m-k-1^{ b};   F_m-2^{ b}   F_m-3^{ b}                         ⋯;            ⋮            ⋮            ⋮            ⋱            ⋮            ⋮; F_m-k+1^{ b}                                      ⋯;   F_m-k^{ b}                                      ⋯               F_m-2k^{ b} ]]

so that the top row consists of the latest sequence values. When multiplied
by the bits, this just generates the next iterate in the sequence.
The upper-diagonal 1's just serve to shift columns over by one, with
each iteration: that is why it's a shift!

The characteristic polynomial of this matrix is, of course, the polynomial
p_n:

    [B-xI]=(-1)^kp_n(x)

Thus, we can trivially conclude that the eigenvalues of B are given
by the roots of p_n(x). This matrix is in lower-Hessenberg
form; this makes it obvious that it's a shift; a finite shift, in
this case.



 §.§ Equivalent labels for orbits


At this point, it should be clear that there are several equivalent
ways of labeling the expressions under consideration. These are recapped
here. Proofs are omitted; they are straight-forward. 



  §.§.§ Orbits


For every given 1<β<2 there is a unique orbit of midpoints
{ m_p} given by m_p=T_β(m_p-1)=T_β^p(m_0)
and m_0=β/2. The orbits are in one-to-one correspondence
with β. The midpoints are the same as the Renyi-Parry seqence;
namely T_β^p(β/2)=t_β^p(1),
recalling here the notation of eqn <ref> and <ref>.



  §.§.§ Orbit encoding


The midpoint generates a unique sequence of bits { b_0,b_1,,⋯,b_k,⋯}
given by the left-right move of the mid-point, as it is iterated.
That is, b_k=Θ(m_k-1/2) so that b_k is
one if the midpoint is greater than half, else b_k is zero. Each
bit-sequence is in one-to-one correspondence with β.



  §.§.§ Monotonicity


The compressor function w(β)=∑_kb_k2^-k
is a monotonically increasing function of β, so that values
of w(β) are in one-to-one correspondence with β. 



  §.§.§ Polynomial numbering


If the orbit is periodic, then there exists a polynomial p_n(z)=z^k+1-b_0z^k-b_1z^k-1-⋯-b_k-1z-1
with k=1+⌊log_2n⌋ being the length
of the orbit. The positive real root of p_n(z) is
β. The integer n is in one-to-one correspondence with the
bit sequence, and with the value of β.



  §.§.§ Integer sequences


If the orbit is periodic, then there exists a sequence of integers
F^{b}, the beta-Fibonacci sequence, that is in one-to-one correspondence
with the finite bit sequence { b} =b_0,b_1,⋯,b_k,
and with the value of β.



  §.§.§ Shift matrix


If the orbit is periodic, then the finite bit sequence { b} =b_0,b_1,⋯,b_k
defines a lower-Hessenberg “golden shift” matrix B, as shown
in eqn <ref>. 



  §.§.§ Summary


To summarize: any one of these: the integer n, the polynomial p_n(x),
the integer sequence F_m^{ b}, the orbit of midpoints
m_p=T^p(β/2), the orbit encoding { b},
the shift matrix B, the value of the compressor function w(β)
and, of course, β itself can each be used as a stand-in for
the others. Specifying one determines the others; they all uniquely
map to one-another. They are all equivalent labels. Fashionably abusing
notation, n≡ p_n(x)≡{ b}≡ F_m^{ b}≡ m_p≡ w(β)≡β≡ B.

An explicit expression relating the orbit encoding and the orbit can
be read off directly from eqn <ref>. Plugging in,

    m_p=T_β^p+1(β/2)=β/2[β^p+1-∑_j=0^pb_jβ^p-j]

for p<k the length of the bit sequence, and m_k=T_β^k+1(β/2)=β p_n(β)/2=0
terminating, since β is the positive root of p_n(x). 

Four of the correspondences given above ask for periodic orbits. Three
of these can be extended to non-periodic orbits in an unambiguous
and uncontroversial way. The extensions are covered in the next two
sections. The fourth is the numbering n of the finite orbits. These
are countable; there is no way to extend the counting number n
to the non-periodic orbits. Ineed, there ae too many: the non-periodic
orbits are uncountable.



 §.§ Infinite-nacci integer sequences 


The beta-Fibonacci integer sequence can be extended to arbitrary (viz.
infinite) length bit sequences, as 

    F_m^{b}=∑_j=1^mb_j-1F_m-j^{b}

starting with F_0^{ b}=1. The sum is always finite,
but one cannot perform it without first knowing at least the first
m bits of the (now infinite) bit-sequence { b}.
The integer sequence still has the desirable property it had before:

    β=lim_m→∞F_m^{b}/F_m-1^{b}

Here, the β value is the one associated to { b}.
So, as before, the real number β and the bit sequence { b}
label exactly the same orbit.

Remarkably, one can be sloppy in how one deals with periodic orbits
with this extention. One has two choices that are equivalent: One
choice is to truncate, so that the bit-sequence ends with all-zeros,
effectively rendering it of finite length. The alternative is to allow
it to continue periodically, forever. Either form results in the same
β-Fibonacci sequence! 

As an example, consider β=1.6, which is close to the golden
ratio, but not quite. It has an infinite non-periodic (non-recurring)
bit-sequence { b} =101010010100101000000100⋯.
The generated integer sequence is F_m^{b}=1,1,1,2,3,5,8,12,20,32,51,82,130,209,335,535,⋯
which undershoots the Fibonacci sequence (12 appears, where we expected
13, and 20 instead of 21, and so on). The ratio of the last two is
535/335=1.597⋯ and the previous is 335/209=1.603⋯
and the ratio of successive elements eventually converges to 1.6.
By comparison, the Fibonacci sequence is generated by the bit-string
1010101010... of alternating ones and zeros. 

The β-Fibonacci representation of the orbits has the remarkable
property that one does not need some a priori mechanism to
know if some orbit is periodic or not. This dual-representation of
periodic orbits is reminiscent of a property commonly seen in Cantor
space 2^ω representations of the real number line, where
the dyadic rationals (which are countable, of course) map to two distinct
bit-sequences (one ending in all-ones, the other ending in all-zeros).
A more general setting for this is given in symbolic dynamics, where
the totally disconnected Bernoulli scheme N^ω can be used
to represent elements of certain countable sets two different ways.
For N=10, one famously has that 1.000...=0.999... as an example.
So likewise here, one can toggle between finite and infinite-periodic
strings. So, given a finite string { b} =b_0,b_1,⋯,b_k-1,b_k
which has, by definition, b_k=1, create a new finite string that
is twice as long: { b^'} =b_0,b_1,⋯,b_k-1,0,b_0,b_1,⋯,b_k
which necessarily has exactly the same beta-Fiboanacci sequence. That
is, F_m^{b^'}=F_m^{b}. Once can repeat this
process ad infinitum, obtaining an infinite periodic string.
The difference between these two is simply the difference between
a less-than-sign, and a less-than-or-equal sign used in the generation
of the orbit, as noted at the very begining of this chapter. We have
proven: finite orbits are exactly the same as infinite periodic orbits,
at least when represented by real numbers and by integer sequences.
Conversely, the difference between using < and ≤ during iteration
is immaterial for describing convergents. 



 §.§ Infinite β-Polynomials


An infinite polynomial is, of course, an analytic function. The goal
here is to extend the finite framework. The definition of the polynomials
above requires a finite bit sequence. This can be extended to an asymptotic
series, by writing first

    p_n(z)=z^k+1(1-b_0z^-1-b_1z^-2-⋯-b_kz^-k-1)

Set ζ=1/z to get

    ζ^k+1p_n(1/ζ)=1-b_0ζ-b_1ζ^2-⋯-b_kζ^k+1

which extends to the holomorphic function

    q^{ b}(ζ)=1-∑_j=0^∞b_jζ^j+1

This is manifestly holomorphic on the unit disk, as each coefficient
is either zero or one. It has a positive real zero, of course: q^{ b}(1/β)=0.
Comparing to eqn <ref>, we see that this is exactly
the same function, or rather, it's negative. Indeed, following the
definition, b_n=d_n(1/2) and so D(β;ζ)=-q^{ b}(ζ). 

This at last provides a foot in the door for correctly describing
the eigenvalues of the β-transfer operator: they are in one-to-one
correspondence with the zeros of q^{ b}(ζ). 



 §.§ β-Hessenberg operator


Extending the golden shift matrix B of eqn <ref>
to an infinite-dimensional operator is a bit trickier. Of course,
one could just declare the matrix elements of the operator to be this-and-such,
but these matrix elements are with respect to what basis? Is the operator
even bounded? The answer to the second question is obviously “no”.

The characteristic equation of B is p_n(β)=0;
the Frobenius-Perron eigenvalue β>1 is too large, although
the k-1=⌊log_2n⌋ other roots are conveniently
arranged near the unit circle, more-or-less equidistant from one another.
The solution is to rescale B by 1/β. The Frobenius-Perron
eigenvalue is now one, and the remaining eigenvalues distributed near
or on a circle of radius 1/β. We may as well take the transpose
as well, so that ℬ_β=B^T/β is in upper-Hessenberg
form. Rescaled in this way, it now seems safe to declare, by fiat,
that the operator ℬ_β is the correct extension
of the matrix B to infinite dimensions. Just to be explicit: given
the bit-sequence { b}, the operatorℬ_β
has the matrix elements

    ⟨ 0|ℬ_β|j⟩    =b_j/β
    ⟨ j+1|ℬ_β|j⟩    =1/β

with all other entries being zero. This is clearly in upper-Hessenberg
form, with the subdiagonal providing the shift.

Comparing to the upper-Hessenberg form of ℒ_β of
eqn <ref>, and the numerical results on it's
eigenvalues, it seems clear that ℬ_β and ℒ_β
must surely be similar. That is, there must be an operator S such
that 


    ℒ_β=S^-1ℬ_βS

The invariant measure ℒ_βρ=ρ is mapped to
σ=Sρ, where ℬ_βσ=σ is the
FP-eigenvector. It is easy to write down σ explicitly: σ=(1,β^-1,β^-2,⋯),
that is, σ_j=β^-j. This is obviously so: the subdiagonal
entries of ℬ_β act as a shift on σ and
the top row is just

    1=∑_j=0^∞⟨ 0|ℬ_β|j⟩σ_j=∑_j=0^∞b_jβ^-j-1=1-q^{ b}(1/β)=1


Although ℬ_β is no more solvable than ℒ_β
in the wavelet basis is, it is certainly much, much easier to work
with. It also re-affirms the ansatz <ref>
for eigenfunctions. To be explicit: if v is a vector satisfying
ℬ_βv=λ v, with vector elements v_j,
then the function

    v(x)=∑_j=0^∞d_j(x)v_j

 is an eigenfunction of the transfoer operator: that is, [ℒ_βv](x)=λ v(x),
or, explicitly:

    1/β[v(x/β)+v(x/β+1/2)]Θ(β/2-x)=λ v(x)

which is just eqn <ref>. So, for λ=1, this
is just v=σ which is just eqn <ref>
for z=1, the invariant measure, as always. But it also says more:
the only solutions to ℬ_βv=λ v are
necessarily of the form v=(1,(λβ)^-1,(λβ)^-2,⋯),
because the subdiagonal forces this shift. To satisfy the the top
row of ℬ_β, one must have that 

    λ=∑_j=0^∞⟨ 0|ℬ_β|j⟩ v_j=1/β∑_j=0^∞b_j/(λβ)^j=λ(1-q^{ b}(1/λβ))=λ

and so the eigenvalue λ is exactly the eigenvalue that solves
the β-series q^{ b}(1/λβ)=0.

This falls short of being a full proof; this line of argumentation
only affirms the ansatz <ref>. To recap:
periodic orbits have an associated shift matrix B; this extends
naturally to a shift operator ℬ_β for non-periodic
orbits. The shift operator has a sufficiently simple form that it's
eigenvectors can be explicitly written down in closed form; they are
necessarily coherent states. The top row of the shift operator defines
a holomorphic function q^{ b} whose zeros correspond
to eigenstates of the shift operator. The holomorphic function is
determined by the binary digit sequence { b}. The
binary digit sequence is obtained from the iterated midpoint, as b_j=d_j(1/2)
where d_j(x)=1 if x<T^n(β/2).
This is enough to prove eqn <ref> holds for
the special value x=0 (for any eigenvalue λ); it
is no enough to show that it holds for any x. For the full proof,
one needs the “fundamental theorem of analytic ergodics”, which
I don't have here.



 §.§ Example eigenfunctions


The above provides a spur to further examine the spectrum of the transfer
operator. Some examples are worked through here.

To recap: eigenstates of the transfer operator correspond with the
zeros of q^{ b}(ζ), or, more precisely,
the zeros for which |ζ|≤1. The reason for this
limitation is that the eigenstates are explictly given by 

    v(x)=∑_m=0^∞d_m(x)ζ^m

for ζ=1/βλ; this is absolutely convergent only for
|ζ|<1. One might hope to analytically extend this
to the entire complex plane, but the extension depends on the digit
sequence d_m(x). We're lacking in a tractable mechanism
to perform this extension.



  §.§.§ Case n=1


Consider first β=φ=1.6180⋯ the golden ratio. The
corresponding finite beta-polynomial is q^{ 11}(ζ)=1-ζ-ζ^2;
the infinite series is 

    q^{ 1010101⋯}(ζ)=1-ζ-ζ^3-ζ^5-⋯=(1-ζ-ζ^2)/(1-ζ^2)

which has a positive real zero at ζ=1/φ and poles at
ζ=±1. The zero corresponds to the FP eigenvalue of one. The
invariant measure is 

    v(x)=∑_m=0^∞d_m(x)/φ^m=φ   0≤ x<1/2
    
    1    1/2≤ x<φ
    
    0    φ≤1

There is a negative real zero at ζ=-φ, but the eigenfunction
summation is not convergent here.



  §.§.§ Case n=2


The n=2 case has the finite bitstring { b} =101
and the periodic bitstring { b} =1001001⋯. The
corresponding finite beta-polynomial is q^{ 101}(ζ)=1-ζ-ζ^3;
the infinite series is 

    q^{ 1001⋯}(ζ)=1-ζ-ζ^4-ζ^7-⋯=(1-ζ-ζ^3)/(1-ζ^3)

which has a positive real zero at ζ=1/β=0.6823⋯ and
three poles on the unit circle. The FP eigenvalue provides β=1.4655⋯.
The invariant measure is 

    v(x)=∑_m=0^∞d_m(x)/β^m=β/β-1   0≤ x<T(β/2)
    1/β-1   T(β/2)≤ x<1/2
    1/β/β-1   1/2≤ x<β
    
    0    β≤1

There are many equivalent ways to write the invariant measure; the
above just selected some representatives from the coset of equivalent
expressions. For example, the third entry could be written as β=1/β(β-1).



  §.§.§ Case n=3


The n=3 case has the finite bitstring { b} =111
and the periodic bitstring { b} =1101101⋯. The
corresponding finite beta-polynomial is q^{ 111}(ζ)=1-ζ-ζ^2-ζ^3;
the infinite series is 

    q^{ 110110⋯}(ζ)=1-ζ-ζ^2-ζ^4-⋯=(1-ζ-ζ^2-ζ^3)/(1-ζ^3)

which has a positive real zero at ζ=1/β=0.5436⋯ and
three poles on the unit circle. The FP eigenvalue gives β=1.8392⋯.
The invariant measure is 

    v(x)=β/β-1   0≤ x<1/2
    β   1/2≤ x<T(β/2)
    1/β-1   T(β/2)≤ x<β
    
    0    β≤1





  §.§.§ Case n=4,6,7


The pattern gets repetitive. There is no case n=5, as this is not
one of the allowed orbits. The bitstrings are those previously listed
in tables; they are { b} =1001 { b} =1101
and { b} =1111. The infinite series is q^{ b⋯}(ζ)=q^{ b}(ζ)/(1-ζ^4).
The zeros are as previously listed. The n=4 plateuas are at 1/β-1[[    β,    1,  1/β, 1/β^2 ]]. The n=6 plateaus are at [[            β/β-1, β^2/(β^2-1)(β-1),                β,    β/(β^2-1)(β-1) ]]. The n=7 plateaus are at [[      β/β-1,          β, β+1/β(β-1),       1/β-1 ]]. Again, the values at the plateaus can be written in many different
ways, given the finite polynomail.



  §.§.§ Case n=16


The n=16 polynomial is the first one to have complex zeros inside
the unit disk. The finite bitstring is { b} =100001
and so the polynomial is q^{ 100001}(ζ)=1-ζ-ζ^6.
The positive real root is ζ=0.7780895986786⋯ and so β=1/ζ=1.28519903324535⋯.
The complex zeros are located at ζ=0.965709509⋯exp± iπ0.2740452363⋯
which corresponds to eigenvalues are λ=0.525107⋯± i0.611100⋯=0.805718⋯exp± iπ0.274045⋯.
The correspnding eigenfunction is shown immediately below.



    < g r a p h i c s >




The order of q^{ b} is six, and this has six almost-plateaus;
they are not quite flat, although they are close to it, presumably
because ζ is close to one.



  §.§.§ The general case


Generalizing from the above, one finds the following:


  * For a period-k orbit, the infinite series is q^{ b⋯}(ζ)=q^{ b}(ζ)/(1-ζ^k). 

  * The first label n for which q^{ b}(ζ)
has a complex zero within the disk is n=16. As a general rule,
it seems that complex zeros inside the disk only appear for β<φ
(I beleive; have not carefully checked. This seems reasonable, as
later chapters show that the region of β<φ behaves very
differently from larger values.)

  * The invariant measure has k plateaus. The plateau boundaries are
given by T^m(1/2) for m={ 0,⋯,k-1}
(so that T^0(1/2)=1/2 and T^1(1/2)=β/2,
and so on).

  * The left-most plateau (of the invariant measure) is at β/(β-1)=∑_n=0^∞1/β^n.

  * The other plateaus appear to be at simple rational functions of β,
but a precise expression is elusive.

To solve the last issue, perhaps one can find tools in Galois theory.
Let ℝ[ζ] be the ring of polynomials in
ζ and consider the quotient ring L=ℝ[ζ]/q^{ b}(ζ).
This L is a field extension of ℝ and so one expects
a Galois group (L/ℝ). The plateaus
of the invariant measure are presumably associated with the group
elements. This seems like a promising direction to go in: perhaps
this is just enough to explain the length of an orbit, the sequence
of points in the orbit, the reason that some polynomials are forbidden
(they don't generate prime ideals), the appareance of Moreau's necklace-counting
function, etc. This remains an unfinished exercise.



 §.§ Factorization


The polynomials factorize. Let r_n denote the real positive root
of p_n(x) – that is, p_n(r_n)=0.
Then one has the factorizations (dropping the subscript on r for
readability)

    p_1(x)=x^2-x-1=(x-r)(x+r-1)=(x-r)(x+p_0(r))

where p_0(x)=x-1. Likewise, there are two order-3
polynomials. They factor as

    p_2(x)=x^3-x-1=(x-r)(x^2+xp_0(r)+rp_0(r))

while


    p_3(x)=x^3-x^2-x-1=(x-r)(x^2+xp_0(r)+p_1(r))

Continuing in this way, there are three order-4 polynomials. They
factor as

    p_7(x)=x^4-x^3-x^2-x-1=(x-r)(x^3+x^2p_0(r)+xp_1(r)+p_3(r))

and

    p_6(x)=x^4-x^3-x^2-1=(x-r)(x^3+x^2p_0(r)+xp_1(r)+rp_1(r))

and (noting that there is no p_5 that occurs in the series)

    p_4(x)=x^4-x^3-1=(x-r)(x^3+x^2p_0(r)+xrp_0(r)+r^2p_0(r))


There's clearly a progression, but perhaps a bit difficult to grasp.
It can be more clearly seen by writing p_n=q_2n+1 and then
writing out 2n+1 in binary. So, once again, from the top:

    p_1(x)=q_11(x)=(x-r)(x+q_1)

where q_1=q_1(r) which adopts the shorthand that
the q polynomials on the right-hand side always have r as an
argument, which can be dropped for clarity. Note also that q_0(r)=r
was already previously observed, in an earlier section. That is, using
the dropped-r convention, q_0=r. Next

    p_2(x)=q_101(x)=(x-r)(x^2+xq_1+q_01)

where, by definition, q_01(x)≡ rq_1(x).
Next,

    p_3(x)=q_111(x)=(x-r)(x^2+xq_1+q_11)

is the second factorization of order 3. For order 4, one has

    p_4(x)=q_1001(x)=(x-r)(x^3+x^2q_1+xq_01+q_001)

where, this time, q_001(x)=xq_01(x)=x^2q_1(x).
Continuing,

    p_6(x)=q_1101(x)=(x-r)(x^3+x^2q_1+xq_11+q_011)

where, by defintion, q_011(x)≡ xq_11(x).
Finally,

    p_7(x)=q_1111(x)=(x-r)(x^3+x^2q_1+xq_11+q_111)

It is worth doing one more, just to clinch that the reversal of the
bit sequence is indeed correct. For this purpose, p_12=q_11001
should serve well. One has

    p_12(x)=q_11001(x)=    (x-r)(x^4+x^3p_0(r)+x^2p_1(r)+xrp_1(r)+r^2p_1(r))
    
    =    (x-r)(x^4+x^3q_1+x^2q_11+xq_011+q_0011)


The general pattern should now be clear. Given one of the admissible
bit sequences b_0b_1b_2⋯ b_k-1b_k and recalling
that b_k=1 always, (and that b_0=1 always) one has 

    p_n(z)=q_b_0b_1b_2⋯ b_k-1b_k(z)=z^k+1-b_0z^k-b_1z^k-1-⋯-b_k-1z-1

which has the factorization, with bits reversed:

    q_b_0b_1b_2⋯ b_k-1b_k(z)=(z-r)(z^k+z^k-1q_b_0+z^k-2q_b_1b_0+z^k-3q_b_2b_1b_0+⋯+q_b_k-1b_k-2⋯ b_1b_0)

where, as already noted, each q is a polynomial in the root r.
Although, notationally, the root r was taken as the real root,
the above factorization works for any root.

The trick can be repeated. Although at first it might seem daunting,
the pattern is uniform: every power of z occurred in the above.
Let s r be some other root. Then

    q_b_0b_1b_2⋯ b_k-1b_k(z)=(z-r)(z-s)(z^k-1+(s+q_b_0)z^k-2+(s^2+sq_b_0+q_b_ob_1)z^k-3+⋯)

The coefficient of the next term being s^3+s^2q_b_0+sq_b_ob_1+q_b_0b_1b_2
and so on. From this point one, this becomes just an ordinary factorization
of polynomials... well, but so was the first step, as well. What made
the first step interesting was that, because the coefficients at that
step were explicitly either zero or one, the corresponding reversal
of the bit sequence became manifest.

One may as well bring this detour to a close. There's nothing particularly
magic in the above factorization, other than the combinatorial re-arrangement
of the polynomial labels. A generic polynomial factorization looks
like the below, for comparison. If 

    p(x)=x^n+1+c_0x^n+c_1x^n-1+⋯ c_n

and if r is a root of p(x) viz p(r)=0
then

    p(x)=    (x-r)(x^n+(r+c_0)x^n-1+(r^2+c_0r+c_1)x^n-2+⋯)
    
    =    (x-r)(x^n+a_0x^n-1+a_1x^n-2+⋯)

with 

    a_k=r^k+1+∑_j=0^kc_jr^k-j

There are some notable values occurring in the factorization. These
are shown in the table below:



o     n     bin     root r     q polynomial     OEIS     root of2     1     11     φ=1.618·      q_1=0.618⋯          4*3     2*2     2*101     2*1.465571·     q_1=0.46557123187676⋯     A088559     q^3+2q^2+q-15-7 
                    q_01=0.68232780382801⋯     A263719      q^3+q-12-7 
     2*3     2*111     2*1.839286·     q_1=0.83928675521416⋯          q^3-2q^2-25-7 
                    q_11=0.54368901269207⋯     A192918     q^3+q^2+q-19*4     3*4     3*1001     3*1.380277·     q_1=r-1          5-7 
                    q_01=0.52488859865640⋯     A072223     q^4-2q^2-q+15-7 
                    q_001=0.72449195900051⋯          2-7 
     3*6     3*1101     3*1.7548776·     q_1=r-1     A075778     q^3+q^2-15-7 
                    q_11=0.32471795724474⋯     silver - 1     5-7 
                    q_011=0.56984029099805⋯          2-7 
     3*7     3*1111     3*1.9275619·     q_1=r-1          5-7 
                    q_11=0.78793319384471⋯          5-7 
                    q_111=0.51879006367588⋯          



As may be seen, some of these constants are already notable for various
reasons. Many are also the real roots of yet other polynomials, of
a not entirely obvious form. (Well, the q_1 polynomials will
always be obvious expansions in binomial coefficients). The suggestion
here is that these are all in turn part of some filigreed partially-ordered
set of intertwining polynomials. Exactly how to express that intertwining
in any sort of elegant or insightful way is not obvious. 



§ ISLANDS OF STABILITY AS ARNOLD TONGUES


The trouble-spots, the eventual fixed-points of the map, can be placed
in one-to-one accordance with the “islands of stability” seen
in the iterated logistic map. They are, in essence, locations where
periodic orbits “could be pasted”, or where they “naturally
would appear”, if the map supported periodic attractors. That is,
the beta shift only supports a single attractor, of period-one at
x=0; there is no “room” for anything more. This is analogous,
in a way, to the phase locked loop, at zero coupling constant. At
finite coupling strength, these “trouble spots” expand out as
Arnold tongues, to have a finite size, visible on the Feigenbaum diagram
for the logistic map as regions where period-doubling is occurring.

The idea here can be illustrated explicitly. Basically, take the natural
sawtooth shape of the map, widen the middle, and insert a slanting
downward line, to create a zig-zag. That is, connect the two endpoints
in the middle of the beta shift, “widening” it so that it has
a finite, not infinite slope, thereby converting the iterated function
from a discontinuous to a continuous one. This can be constructed
directly: given some “small”, real ε>0, define the
piecewise-linear ε-generalization of the map <ref>
as 

    T_β,ε(x)=β x    0≤ x<1/2-ε
    β/4-β(1/4-ε)w    1/2-ε≤ x<1/2+ε
    β(x-1/2)    1/2+ε≤ x≤1

where w is just a handy notation for a downward sloping line:

    w=2x-1/2ε

Observe that w=1 when x=1/2-ε and that w=-1
when x=1/2+ε so that w just smoothly interpolates
between +1 and -1 over the middle interval. The additional factors
of β/4-β(1/4-ε)w just
serves to insert the downward slope smack into the middle, so that
the endpoints join up. The results is the zig-zag map, illustrated
in the figure below



    < g r a p h i c s >




In the limit of ϵ→0, one regains the earlier beta shift:
lim_ε→0T_β,ε=T_β, as the slope
of the middle bit becomes infinite. The middle segment is a straight
line; it introduces another folding segment into the map. This segment
introduces a critical point only when ε is sufficiently
large, and β is sufficiently small, so that its slope is less
than 45 degrees (is greater than -1). When this occurs, a fixed point
appears at x=1/2. A sequence of images for finite ε
are shown in figure <ref>.



The appearance of islands of stability in the Feigenbaum attractor
is due to the presence of a fixed point at any parameter value. In
order to “surgically add” islands of stability to the beta transform,
the middle segment interpolation must also have a critical point at
“any” value of ε. To achieve this, consider the curve


    D_β,ε(x)=β x    0≤ x<1/2-ε
    β/4-β(1/4-ε)g(w)    1/2-ε≤ x<1/2+ε
    β(x-1/2)    1/2+ε≤ x≤1

where the straight line has been replaced by a soft shoulder 

    g(w)=1-2cosπ/4(1+w)

and w is the same as before. This is scaled so that its a drop-in
replacement for the straight line: g(1/2-ε)=1
and g(1/2+ε)=-1. A cosine was used
to create this soft shoulder, but a parabola would have done just
as well. It is illustrated above, with the label “soft map”.

This map also interpolates between the left and right arms of the
beta transform, forming a single continues curve. The curve is smooth
and rounded near 1/2-ε x, having a slope
of zero as x approaches 1/2-ε from above.
This introduces a critical point near 1/2-ε.
Notice that there is a hard corner at 1/2+ε.
The interpolation is NOT an S-curve! A sequence of images for finite
ε are shown in figure <ref>.



Two more variant maps can be considered. Both replace the center piece
with symmetrical sinuous S-shaped curves, but in different ways. Consider

    S_β,ε,σ(x)=β x    0≤ x<1/2-ε
    β/4-σβ(1/4-ε)sinπ/2w    1/2-ε≤ x<1/2+ε
    β(x-1/2)    1/2+ε≤ x≤1

and

    H_β,ε,p,σ(x)=β x    0≤ x<1/2-ε
    β/4-σβ(1/4-ε)(x-1/2)|w|^p   1/2-ε≤ x<1/2+ε
    β(x-1/2)    1/2+ε≤ x≤1

The S_β,ε(x) replaces the central segment with
a softly-rounded segment, containing two critical points: near 1/2-ε
and near 1/2+ε, where the curve flattens out
to a zero slope. When σ=+1, the map as a whole is continuous.
When σ=-1, the map consists of three discontinuous pieces.
Different values are explored in figure <ref>.



The H_β,ε,p,σ(x) replaces the central segment
with a segment that has a kink in the middle, when p>1. Note that
H_β,ε,1,1(x)=T_β,ε(x). Here, x
is the sign of x. The general shape of H_β,ε,p,σ(x)
is shown above, labeled as the “kink map”. The location of the
kink in H is always centered; an off-center kink, as depicted in
the figure, is explored below. The bifurcation diagrams for H are
illustrated in figure <ref>.



To summarize: the “trouble spots” don't “just” break the ability
to create a Hessenberg basis at certain values of β: they are
more “fundamental” than that: they indicate the regions where
(“phase-locked”) periodic orbits can be made to appear. The last
sequence of images, shown in figure <ref>
indicate that the islands of stability need NOT consist of the period-doubling
sequences seen in the Feigenbaum map. This is made explicit in figure
<ref>, which shows a zoom by a factor of thirty.



Another interesting visualization is a Poincaré recurrence plot. The
islands of stability should manifest as Arnold tongues<cit.>.
These are shown in figures <ref> and
<ref>.





To intuitively understand the location of the islands (the location
of the Arnold tongues), its easiest to examine a map with a kink in
it, whose location is adjustable. 


    H_β,ε,α,σ(x)=β x    0≤ x<1/2-ε
    β/4-σβ(1/4-ε)h_α,p   1/2-ε≤ x<1/2+ε
    β(x-1/2)    1/2+ε≤ x≤1

with 

    h_α,p(x)=α+(1-α)|w|^p   x<1/2
    α-(1+α)|w|^p   1/2≤ x

As before, h_α,p(x) is designed to interpolate
appropriately, so that h_α,p(1/2-ε)=1
and h_α,p(1/2+ε)=-1. The location
of the kink is now adjustable: h_α,p(1/2)=α.
Iterating on this map results in figures that are generically similar
to those of figure <ref>, except that
this time, the location of the islands is controllable by the parameter
α. Roughly, to first order, the primary series of islands
are located at √(2/(1-α)); as before, these
islands do not allow period-doubling to take place. 

To get islands with period doubling, one needs to re-create the “soft
shoulder” of eqn <ref>, but at a variable location.

Thus, the above presents a general surgical technique for controlling
both the general form of the chaotic regions, the location of the
islands of stability, and what appears within the islands.

Note to reader: I suspect that the above observations have been previously
discovered by others, and might even be “well known”, viz.
presented in some pop-sci literature on fractals. However, I am not
aware of any references discussing this topic. If you, dear reader,
know of such references, please drop me a line at the posted email
address.

Exercise left to the reader: the above arguments should be sufficient
to fully demonstrate that the circle map, which is well-known to exhibit
phase locking regions called Arnold tongues, is topologically conjugate
to the fattened beta shift T_β,ε. Or something
like that. In a certain sense, this can be argued to be a “complete”
solution, via topological conjugacy, of the tent map, the logistic
map and the circle map. This is a worthwhile exercise to actually
perform, i.e. to give explicit expressions mapping the various regions,
as appropriate. 

Essentially, the claim is straight-forward: topologically, all chaotic
parts of a map correspond to folding (as per Milnor, 1980's on kneading
maps), into which one may surgically insert regions that have cycles
of finite length. The surgical insertion can occur only at the discontinuities
of the kneading map. It almost sounds trivial, expressed this way;
but the algebraic articulation of the idea would be worthwhile.



§ MISCELLANEOUS UNFINISHED IDEAS


An ad-hoc collection of half-finished thoughts.



 §.§ Multiplicative Shifts


A multiplicative shift is a shift assembled as an product of functions.
The most famous of these is the generating function for integer partitions

    P(z)=∏_n=1^∞1/(1-z^n)

Similarly products occur for the necklace counting functions, most
famously the cyclotomic identity

    1/1-β z=∏_j=1^∞(1/1-z^j)^M(β,j)

where M(β,j) the necklace polynomial.

A far more obscure product expresses the Minkowski measure<cit.>,
given as

    ?^'(x)=∏_n=0^∞A^'∘ A^n(x)/2

with 

    A(x)=x/1-x   0≤ x<1/2
    2x-1/x   1/2≤ x≤1

with A^' being the derivative of A and A^n being
the n'th iterate. The Minkowski measure integrates to the Minkowski
Question Mark function ?(y)=∫_0^y?^'(x)dx;
it is the prototypical “multi-fractal measure” (although there
really is nothing “multi-” about it; the “multi-” prefix stems
from a misunderstanding of its multiplicative invariance). The product
structure indicates that the Minkowski measure is a Gibbs measure,
viz arising from an invariant Hamiltonian on a one-dimensional lattice.

The figure <ref> suggests that a similar
product can be constructed from the midpoint sequence, namely 

    ∏_p=0^∞4m_p(β)/β

for the midpoints m_p(β)=T_β^p(m_0).




 §.§ Midpoints, revisited


The midpoints are defined above as m_0=β/2, so that m_p=T_β(m_p-1)=T_β^p(m_0)
with T_β(y) the beta shift map of eqn <ref>.
Almost all literature uses the beta transform t_β(x)
of eqn <ref> instead. The midpoint sequence and
the iterate t_β^p(1) are closely related:

    2m_p1=t_β^p+1(1)

Although related, they are not the same. The difference is a sequence
of bits:

    c_p=2m_p-t_β^p+1(1)

Note that c_p∈{ 0,1} always. Note that

    β=∑_p=0^∞c_p/β^p

which is not entirely obvious!



 §.§ Rauzy Fractals


Given a polynomial, one has an associated finite matrix, in Hessenberg
form, that, iterated upon, generates a sequence. The projection of
that sequence to a non-expanding orthogonal plane is a Rauzy fractal.
What are the corresponding Rauzy fractals for this situation?

How about the general iterated sequence (e.g. the sequence of midpoints)?
Is this space-filling, or not? 



§ BERGMAN (HESSENBERG) POLYNOMIALS


Given a matrix operator in Hessenberg form, it can be interpreted
as a right-shift on the space of polynomials. Such polynomials form
an orthogonal basis for a certain kind of Hilbert space, called the
Bergman space. They are studied in applied mathematics, as they are
orthogonal over some measure on the complex plane. The Hessenberg
operator is a generalization of the better-known case of the Jacobi
operator, which has it's own extensive theory, including spectra and
scattering, and is important for several exactly solvable non-linear
models in physics, including the Toda lattice<cit.>. The
Hessenberg operator presumably has an equally rich theory, but it
does not appear to be currently known; the breadth and scope of existing
publications is limited.

The general framework for the Hessenberg polynomials is sketched below,
including a fast and informal definition of Bergman space (the space
on which the polynomials are orthogonal). The Hessenberg matrix is
explicitly solvable on the left, and can be explicitly brought into
a from that exhibits the right-shift operator. In the general theory,
the change of basis from the shift operator to the Hessenberg matrix
is known to be the Cholesky decomposition of a moment matrix, and
specifically, the moments of the measure on which the polynomials
are orthogonal.

There are two Hessenberg operators in this text: the operator ℒ_β
in the wavelet basis, and the operator ℬ_β generated
from the midpoint orbits. The second is already obviously a shift,
and so everything below follows “trivially” from it. The first
form is is numerically and analytically difficult. Needless to say,
the section below treats the firs rather than the second. XXX TODO
this should be fixed, as ℬ_β is both simpler and
more enlightening overall. Later ... 

Working backwards from the beta shift, the first asymptotic term in
the measure can be extracted. For β>φ, it appears to
be a Dirac delta (point mass) located at z=1 on the complex plane,
with a blancmange-like fractal curve giving the weight. For β<φ,
it appears to be the derivative of the Dirac delta, with a different
blancmange-like fractal curve giving the weight.



 §.§ Bergman Space


Given a matrix operator in Hessenberg form, it can be interpreted
as a right-shift on the space of polynomials. That is, given an unreduced
Hessenberg matrix with matrix entries A_ij, one can write a recurrence
relation that defines a sequence of polynomials as

    zp_n(z)=∑_k=0^n+1A_knp_k(z)

with p_0(z)=1. This relation is easily solvable in
closed form, as the recurrence relation terminates in a finite number
of steps. 

One important property of these polynomials is that the zeros of p_n(z)
correspond to the eigenvalues of the n× n principle submatrix
of A. Numeric exploration of these polynomials confirms the previous
results on eigenvalues obtained from direct diagonalization: the zeros
of the p_n(z) seems to lie mostly near the circle
of radius 1/β, distributed uniformly over all angles.

If all of the sub-diagonal entries obey A_n+1,n>0, then the polynomials
form an orthonormal basis for Bergman space. That is, there exists
a domain in the complex plane on which the polynomials provide a basis
for a Hilbert space of holomorphic functions on that domain<cit.>.
That is, one has the orthogonality relation

    δ_mn=∫_Dp_m(z)p_n(z)dμ(z)

for some domain D⊂ℂ of the complex plane, and some
(Borel) measure dμ on that domain. 

The matrix A can be interpreted as an operator with a continuous
spectrum. To do this, fix a certain, specific value of z=c a constant,
and then notice that p⃗=(p_n(z))_n=0^∞
is a vector having the property that A^Tp⃗=zp⃗. That
is, p⃗ is a left-eigenvector of A; equivalently, a right-eigenvector
of its transpose A^T. Clearly, the spectrum is continuous on
the domain D.

The matrix operator A can also be interpreted as a right-shift
on Bergman space. To do this, define

    𝒜(w,z)=∑_k=0^∞∑_n=0^∞p_k(w)A_knp_n(z)

Then, given some holomorphic function f(z) decomposed
in terms of the polynomials, so that f(z)=∑_na_np_n(z),
one has that

    [𝒜f](w)=    ∫𝒜(w,z)f(z)dμ(z)
    
    =    ∑_k∑_np_k(w)A_kna_n
    
    =     w∑_na_np_n(w)
    
    =     wf(w)

That is, given a sequence (a_0,a_1,a_2,⋯),
the Hessenberg matrix acts as a right-shift, mapping it to the sequence
(0,a_0,a_1,⋯). 

This is perhaps a bit silly, as one could instead just perform the
same manipulation without the f(z), by observing that,
formally,

    𝒜(w,z)=w∑_k=0^∞∑_n=0^∞p_k(w)p_n(z)


The above treatment is breezy and “formal”, paying no heed to
summability, convergence or responding to any questions about what
spaces the various vectors may live in. This is as appropriate, since
the task here is to discover which spaces are the appropriate ones,
when the Hessenberg matrix arises from the beta shift.

Notice that the word “operator” is a bit mis-used, here, as a
vague synonym for “infinite-dimensional matrix”. Properly, the
word “operator” should be reserved for an infinite-dimensional
matrix acting on some given space, having general properties that
are independent of the basis chosen for that space. So far, that might
not be the case here: the infinite-dimensional matrices here might
not be bounded operators; they might not even be continuous, viz.
we have not ruled out the possibility that the space of interest is
some Fréchet space or some more general topological vector space.
It is well known that operators on such spaces can have “unexpected”
discontinuities, unexpected in that they are not seen in ordinary
Banach spaces.

At any rate, if polynomials obtained from the beta shift are orthogonal
on some domain D⊂ℂ that is the support of some measure
dμ, it is not at all clear what this measure might be. They are
certainly not orthogonal on the unit disk, with uniform measure.

Notice also that the above treatment seems to be a special case of
a more general principle: when an operator has a continuous spectrum,
it can sometimes be interpreted as a right-shift. That is, given some
arbitrary operator ℋ, then if one has that ℋf=λ f
and λ takes arbitrary values λ∈ D⊂ℂ,
then ℋ can be taken to be a right-shift operator, provided
that f=f(λ) can be decomposed into a set of orthogonal
polynomials in λ. 



 §.§ Beta Bergman Shift


The primary question for this section is whether the β-transform
transfer operator, in the Hessenberg basis, can be considered to be
a Bergman shift.

To obtain the orthogonal polynomial basis, one must satisfy the constraint
that A_n+1,n>0 for the matrix elements A_kn=⟨ k| L_β|n⟩
of eqn <ref>. Numeric exploration indicates
that this is indeed the case, with the sub-diagonal entries all positive
(none are zero), and all tend to have the same value, with sporadic
exceptions. These are shown in figure <ref>. 



Can one find a domain on the complex plane that would have such Bergman
polynomials? The references<cit.> provide a technique
for doing so, provided that the matrix is asymptotically Toeplitz.
That is, if the diagonals of A_ij have nice limits, that lim_n→∞A_n-k,n
exists for fixed k, then a Jordan arc bounding a domain on the
complex plane can be found. The figure <ref>
indicates that this limit does not exist, in the strict sense: the
values bounce away from an obvious limit point indefinitely. Exactly
what this implies is unclear. Perhaps it is possible to extend the
results of <cit.> to matrices that are where the diagonals
merely have an accumulation point, as opposed to a well-defined limit? 

Based on numeric exploration, it appears that the domain is the unit
disk. That is, A^Tp⃗=zp⃗ holds for |z|≤1.



 §.§ Bergman Alternative


The Bergman polynomials of eqn <ref> define
an orthonormal basis for some region of the complex plane. For the
square-integrable norm, this basis is the basis of a Hilbert space,
and specifically, that of a reproducing kernel Hilbert space. 

Yet, something funny happens on the unit disk. Let p_m(z)
be the polynomials, and for some sequence of coefficients { a_n},
consider a generic function 

    f(z)=∑_k=0^∞a_kp_k(z)

Consider the case where the { a_n} are a right-eigenvector
of the Hessenberg operator, that is, where 


    ∑_m=0^∞A_kma_m=λ a_k

Substituting into the above, one has 

    f(z)=∑_k=0^∞1/λ∑_m=0^∞A_kma_mp_k(z)=z/λ∑_m=0^∞a_mp_m(z)=zf(z)/λ

There are two alternatives to solving this; either f(z)=0
or z=λ. Since this is a reproducing kernel Hilbert space,
then if z=λ is part of the domain of the Bergman space, then
one must conclude that f(z)=0 everywhere. That is, right-eigenvalues
of A correspond to functions f(z) that are vanishing.
To invent a new name, by analogy to the Fredholm alternative, perhaps
this can be called the Bergman alternative.

Numerical exploration indicates that, for the matrix elements of eqn,
<ref>, the function f(z) vanishes
inside the unit disk |z|<1, and is undefined (infinite)
outside of it.



 §.§ Left Factorization


Suppose one is given an (arbitrary) sequence of polynomials (p_n(z))_n=0^∞,
such that the order of p_n is n. Then each individual polynomial
can be expanded as or β>φ,

    p_n(z)=∑_k=0^np_nkz^k

This defines an infinite matrix 𝒫=[p_nk],
provided that the coefficients are extended so that p_nk=0 whenever
k>n. This matrix is manifestly lower-triangular. Writing vectors
z⃗=(z^n)_n=0^∞ and p⃗=(p_n(z))_n=0^∞
as before, the above is just the matrix equation 

    p⃗=𝒫z⃗


Consider now the case where the polynomials were constructed from
some irreducible Hessenberg matrix A. The earlier observation that
A^T is a shift, namely, that A^Tp⃗=zp⃗ can now
be written as

    A^T𝒫z⃗=z𝒫z⃗=𝒫zz⃗=𝒫𝒦z⃗

In the above, the z without the vector notation is just a scalar,
and thus commutes (trivially) with 𝒫. Its eliminated
by explicitly making use of the right-shift (Koopman) operator, which,
in this basis, is

    𝒦=[[ 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 ⋱; 0 0 0 0 ⋱ ]]

Since 𝒫 is lower-triangular, it is invertible on the
right, that is, the inverse 𝒫^-1 exists, and so one
is left with 

    𝒫^-1A^T𝒫=𝒦

The irreducibility of A is important, here; non-zero entries on
the sub-diagonal are required, else trouble ensues.

Rearranging, this provides an explicit decomposition of A into
triangular matrices:

    A^T=𝒫𝒦𝒫^-1

Taking the transpose, this gives

    A=[𝒫^-1]^T𝒦^T𝒫^T

with 𝒫^T and [𝒫^-1]^T both
being upper-triangular, and 𝒦^T being the left-shift.

This system is solvable. Given some matrix A in Hessenberg form,
the matrix elements of 𝒫 can be computed recursively,
in a finite number of steps (i.e. in closed form), directly from <ref>.
The explicit expression is


    A_n+1,np_n+1,j=p_n,j-1-∑_k=0^nA_knp_kj

The starting conditions are p_00=1. To handle the j=0 case
in the above, set p_n,-1=0. 

Because 𝒫 is lower triangular, its inverse 𝒫^-1≡ℛ=[r_kn]
can be obtained explicitly. Along the diagonal, one has r_nn=1/p_nn
while the lower triangular form means r_kn=0 for k<n. For
the remaining entries m<n, one has 

    0=∑_k=m^np_nkr_km

This can be solved in a finite number of iterations on 

    p_nnr_nm=-∑_k=m^n-1p_nkr_km


The above avoids questions of convergence, or any notion of the spaces
on which the matrices or operators might act. The norm to be used
for z⃗ and p⃗ is not specified. This is appropriate
at this stage: it is the algebraic manipulations that are interesting,
at this point, rather than the spaces on which the matrices/operators
might act. One can invent several kinds of norms that might be applicable,
but there is no particular reason to believe that p⃗ might
have a finite norm. Likewise, 𝒫 may not have a finite
norm. For the case of the Hessenberg operator originating with the
beta shift operator, it does not; the individual matrix elements p_nm
increase without bound. That is, 𝒫 is an infinite matrix,
but it is not clear that it is also an operator. If it is, it is certainly
not a compact operator.

Some of the poor behavior can be brought under control by factoring
𝒫=𝒟𝒩 with 𝒩 being unitriangular
(all ones on the diagonal) and 𝒟 a diagonal matrix, with
entries [𝒟]_nk=p_nnδ_nk. With this
factorization, one may then write

    𝒩^-1A^T𝒩=𝒟𝒦𝒟^-1

so that 𝒟𝒦𝒟^-1 has off-diagonal matrix entries [𝒟𝒦𝒟^-1]_nk=δ_n+1,kp_nn/p_kk.
This is a rescaling of the shift [𝒦]_nk=δ_n+1,k.
The scaling factor is exactly the sub-diagonal of the Hessenberg.
That is, p_nn/p_n+1,n+1=A_n+1,n. The polynomials 𝒩z⃗
are monic.



 §.§ Beta-transform factoids


An assortment of observations follow, for the case of the beta shift.

First, the matrix entries of 𝒫 grow in an unbounded fashion.
It appears that p_nn∼𝒪(β^n); the
ratio p_nn/β^n is depicted in figure <ref>.



Experimentation reveals two different regimes of behavior, depending
on whether or not β<φ=(1+√(5))/2 the Golden
ratio. Exactly why there are two different regimes is unclear. Earlier
sections motivated the reason for the appearance of the golden mean;
why this shows up dramatically, as it does here, is unclear (to me).

One such result is that when β<φ, then the sum over columns
of the Bergman operator vanishes. That is,

    ∑_k=0^∞p_nk=δ_n0

This implies that every polynomial p_n(z) has a zero
at z=1 (except for p_0(z)=1) when β<φ. 



 §.§ Decaying Eigenfunctions


The matrix mechanics developed in the previous sections can be used
to perform asymptotic expansions that rapidly converge to decaying
eigenfunctions. This works most simply for the case of φ<β.
TODO Write these down. TODO flesh out. Basically, write a vector w⃗
with elements w_n=ω^n for 1<|ω| so
that this is divergent. Then write the formal vector a⃗=[𝒫^T]^-1w⃗
which is formally divergent, but can be truncated in finite dimensions,
and renormalized to be of unit length. Doing so provides an eigenfunction
of A. The associated eigenvalue is 1 when β<φ but
is less than 1 when φ<β (and in fact, the eigenvalue
is exactly that depicted in figure <ref>).
TODO graph some of these, explore more thoroughly, address the issues
of formal divergence.



 §.§ Moment Matrix


When the Hessenberg matrix is derived from measures on the complex
plane, it takes the form of ℳ=ℛℛ^T with ℛ=𝒫^-1,
so that ℛ is the Cholesky decomposition of ℳ.
This matrix is manifestly symmetric: ℳ=ℳ^T.
Direct observation shows that it is almost positive-definite: one
finds that [ℳ]_ij>0 for all i,j except
for [ℳ]_00=0. This result can be strengthened:
when β<φ, then [ℳ]_ij>1 for
all i,j except for [ℳ]_00=0 and [ℳ]_0n=[ℳ]_n0=1.
But, for β>φ, one finds that [ℳ]_00=0
and [ℳ]_01=[ℳ]_10=[ℳ]_11=1,
while all the rest obey 0<[ℳ]_ij<1.

In the standard literature, ℳ is usually obtained from
some moment matrix, viz, for the integral ∫z^mz^ndμ(z)
for some measure dμ(z). Might that be the case here?
Taking the time to numerically characterize the matrix, one finds
that the ratio of successive rows (or columns as its symmetric) very
quickly approaches a limit lim_n→∞,[ℳ]_nm/[ℳ]_n-1,m=C(β)
for some constant C that depends only on β but not on m.
The limit C(β) is graphed in figure <ref>.



For β<φ, it appears that lim_n→∞,[ℳ]_nm=B(β)
a constant, independent of m. This limiting value B(β)
is graphed in figure <ref>.



The asymptotic behavior of the matrix [ℳ]_ij
can be obtained as a moment matrix on point sources. A delta function
located at z=C for real C has the moments

    C_mn=    ∫z^mz^nδ(z-C)dz
    
    =    ∫ r^mr^nδ(r-C)rdr∫δ(θ)e^-imθe^inθdθ
    
    =     C^m+n+1

Thus, for φ<β, the asymptotic behavior of [ℳ]_ij
is given by the distribution A(β)δ(z-C(β)).
What is A(β)? This is graphed in figure <ref>.



What about β<φ? A limiting constant distribution can
be obtained from a derivative point mass located at z=1. That is,

    D_mn=    ∫z^mz^nδ^'(z-1)dz
    
    =    ∫ r^mr^nδ^'(r-1)rdr∫δ(θ)e^-imθe^inθdθ
    
    =     1

so that the asymptotic behavior of [ℳ]_ij
for β<φ is given by the distribution B(β)δ^'(z-1).
The prime superscript here means derivative, viz, in colloquial language,
δ^'(z)=dδ(z)/dz.



§ THE JACOBI OPERATOR


Given a Borel measure on the real number line, one can find a sequence
of polynomials that are orthonormal with respect to that measure.
These polynomials p_n(x) are coupled together by a
three-term recurrence equation

    xp_n(x)=a_n+1p_n+1(x)+b_np_n(x)+a_np_n-1(x)

with p_0(x)=1 and p_-1(x)=0. This
recurrence relation can be taken to be an operator, known as the Jacobi
operator 𝒥, acting on vectors consisting of the polynomials
p(x)={ p_n(x)} so that

    [𝒥p](x)=xp(x)

so that p is an eigenvector of 𝒥 with eigenvalue x.
The two sequences of coefficients { a_n} and { b_n}
form three diagonals of the operator, with { a_n}
running down the center, and { b_n} the two diagonals
on either side<cit.>.

Given that the invariant measure for the β-transform, given
by eqn <ref> and visualized in figure <ref>
is a Borel measure, it seems reasonable to ask: what is the corresponding
Jacobi operator? How can the sequence of polynomials be understood? 

Szegő polynomials w.r.t. dμ are a set of orthogonal polynomials
on the unit circle. Applying a Cayley transform gives the Schur functions,
obeying a rational recurrence relation solvable via continued fractions.
Hmmm.

And then there is Favard's theorem...



 §.§ Moments


Construction of the polynomial sequences require moments. Since the
invariant measures (and all of the eigenfunctions) are linear combinations
of the Hessenberg basis functions, it suffices to compute the moments
for these. Since the basis functions are piece-wise constant, and
have an explicit expression given by eqn <ref>,
the moments can also be given explicit expression:

    ∫_0^1x^n-1ψ_p(x)dx=C_p/n[m_p^n-m_l^n/m_p-m_l-m_u^n-m_p^n/m_u-m_p]

with the midpoint m_p and the lower and upper midpoints m_l<m_p<m_u
defined just as before. Clearly, the moments rapidly get small as
n→∞. Likewise, for fixed n, these also rapidly get small
as p→∞. 



§ THE MULTIPLICATION OPERATOR


The difficulties presented in the previous section suggests that studying
the multiplication operator might be simpler. Multiplication by β
is given by 

    M_β(x)=β x

The corresponding transfer operator is 

    [ℳ_βf](y)=1/βf(y/β)

The multiplication operator, superficially, in itself, is not terribly
interesting; it simply rescales things. It does not generate fractals,
at least, not if one confines oneself to real numbers and the canonical
topology on the real-number line. If instead one works with the product
topology on 2^ω, then the multiplication operator becomes
rather complicated and difficult to analyze. In this sense, it is
promising: it avoids the overt complexity of the logistic map, the
tent map and the beta shift, yet still has a complicated behavior
in the product topology. In particular, the multiplication of two
numbers appear to involve chaotic dynamics of the carry bit.



 §.§ Beta-shift, Revisited


The beta shift of eqn <ref> takes a simple form when
reinterpreted on bit-strings: it is the concatenation of multiplication,
followed by a left-shift. Given a bit-string (b_n)=0.b_0b_1b_2⋯
denote its left-shift by U given by

    U(0.b_0b_1b_2⋯)=0.b_1b_2⋯

which, for real numbers, corresponds to 

    U(x)=
    2x    0≤ x<1/2
    
    2x-1    1/2≤ x≤1

which is none-other than the Bernoulli shift of eqn <ref>
with a change of notation. The beta shift is then

    T_β(x)=M_β(U(x))

so that the iterated beta shift is an alternation between a left-shift
and a multiplication. The act of discarding the most significant bit
(the MSB) with each left-shift alters the dynamics of iterated multiplication.

This suggests that studying multiplication and the multiplication
operator might provide fruitful insight into the beta shift.



 §.§ Monomial Eigenfunctions


Some properties of the multiplication operator can be guessed at directly.
Obviously, f= is a decaying/growing eigenfunction,
depending on whether β>1 or not. That is, one should imagine
f= as a uniform distribution of dust; with each iteration,
it is spread either farther apart (β>1) or bunched closer together
(β<1).

Clearly, f(x)=x^n is an eigenfunction, with eigenvalue
1/β^n+1. If one considers multiplication only to operate
on the positive real-number line, then n need not be an integer.
In other words, the multiplication operator has a continuous spectrum
in this situation.

If the domain of the operator is extended to functions on the non-negative
real-number line, then n must be positive, as otherwise f(0)
is ill-defined. But if n is positive, then (for β<1) the
multiplication operator only has eigenvalues greater than one, which
is not, in general, very desirable.

If the domain of the multiplication operator is extended to the entire
real-number line, then n is forced to be an integer, in order to
avoid issues due to multi-valued functions. Extending the domain to
the complex plane leads us astray, and so we will not go there.



 §.§ A Fractal Eigenfunction


The compressor function is also an eigenfunction. It was previously
observed in eqn <ref> that 

    _β(x/β)=1/2_β(x)

whenever 1<β≤2 and 0≤ x<1 and so, _β
is potentially be an eigenfunction of ℳ_β with
eigenvalue 1/2β, provided that it is extended to arguments
1<x. This can be done as follows. Define the extended function,
valid for 0≤ x<∞ and for 1<β≤2 as

    _β(x)=_β(x)    0≤2x<β
    
    2_β(x/β)    β≤2x<β^2
    
    4_β(x/β^2)    β^2≤2x<β^3
    
    2^n_β(x/β^n)    β^n≤2x<β^n+1

The extension is performed simply by treating the self-similarity
as a recurrence relation, which can be iterated to move the argument
into a region where the original definition was sufficient. In essence,
one applies a right-shift operator to reduce the argument. Since the
multiplication operator is odd about x=0, on can trivially extend
this to negative x by defining _β(-x)=-_β(x). 

Note that the original _β(x) also had
a translation symmetry: the upper half was equal to the lower half.
This translation symmetry has been lost, since after all, multiplication
does not preserve translation.

The ecpr function is not square integrable; it does not have an L_p-norm
for any p; and this is no surprise, as its hard to imagine how
it could be otherwise, for a function to be self-similar under scaling. 



 §.§ A Generic log-periodic Eigenfunction


Inspired by the above, its should be clear how to build a generic
eigenfunction. Let g(x) be some arbitrary function,
defined on the interval 1≤ x<β (given some fixed 1<β).
Define its extension as

    g_w^'(x)=w^ng(x/β^n) β^n≤ x<β^n+1

This has, by construction, the self-similarity relation g_w^'(β x)=wg_w^'(x)
and so is an eigenfunction with eigenvalue w/β:


    [ℳ_βg_w^']=w/βg_w^'

This function is merely log-periodic; its not fractal. Perhaps its
silly to illustrate this; it should be obvious, but just in case its
not, the figure below shows such a function, for β=1.6 and
w=0.8. It is an eigenfunction of ℳ_1.6 with eigenvalue
of 1/2.



    < g r a p h i c s >






There doesn't seem to be anything particularly interesting with this
particular game. There's a simple explanation for this: The multiplication
operator is generating a free monoid in one generator (the iteration
itself), whereas fractals require at least two generators of self-symmetry.
The (usually) free interaction of multiple generators is what forces
the fractal to appear.



Note that the _β function constructed above is
a special case of this: It's self-similar, but the property that made
it interesting, as a fractal, was erased in the construction. As before,
note that g_w^'(x^n) is an eigenfunction
with eigenvalue 1/β w^n (for integer n).




 §.§ Haar Basis Matrix Elements


The Haar basis matrix elements for the beta shift proved to be a bit
unwieldy and not terribly useful. The corresponding matrix elements
for the multiplication operator have the same general essence, but
are slightly simpler and shorter to write down. In all other respects,
they still have the same tractability issues. 

The multiplication operator ℳ_β has matrix elements
in the standard Haar basis:

    ⟨ mi| M_β|nj⟩ =    ∫_-∞^∞h_mi(x)[ℳ_βh_nj](x)dx
    
    =    2^(m+n)/2/β∫_-∞^∞h(2^mx-i)h(2^nx/β-j)dx

Instead of confining oneself to the unit interval, here it is convenient
to consider the entire real-number line, and thus that is the range
of the integral. Likewise, i and j an be any integers, positive
or negative. As before, matrix elements vanish unless 

    [i/2^m,i+1/2^m]∩[β j/2^n,β(j+1)/2^n]∅

This holds in three cases: where one of the intervals contains an
edge transition (left, middle or right) of the other interval, without
also containing the other two.



 §.§ The Shift and Add algorithm


One can model the multiplication of real numbers with a number of
different algorithms applied to bit strings. One of the simplest such
algorithms is the shift-and-add algorithm, described here. Its just
elementary-school long-form multiplication, applied to the binary
expansions of the numbers. 

There's a point worth laboring on: a bit string representing a real
number is not the same thing as the real number. There are more bit-strings
than there are real numbers. Most famously, the two bit strings 0.0111⋯
and 0.1000⋯ are two obviously distinct bit-strings, but they
represent the same real number: one-half. All real numbers of the
form j/2^n (the “dyadic rationals”) will always have dual
representations; all other real numbers have a single, unique representation.
These correspond to the “gaps” in the Cantor set, or, equivalently,
neighboring infinite branches in the finite binary tree. Bit-strings
are not real numbers. They're just a usable model of them. The usability
is somewhat limited; its OK for working with individual points, but
fails miserably for the topologies: the canonical topology on the
reals is sharply different than the product topology on 2^ω.

The goal is to compute the product Kx with 0≤ K≤1 and 0≤ x≤1
so that the product is 0≤ Kx≤1. Both K and x are represented
by their binary expansions. Let the binary expansions be

    x=0.b_0b_1b_2⋯=∑_n=0^∞b_n2^-n-1

and


    K=0.c_0c_1c_2⋯=∑_n=0^∞c_n2^-n-1

where the b_n and c_n are either 0 or 1, always.

Define s_0=0 and s_n+1 to be the non-negative integer

    s_n+1=b_nc_0+b_n-1c_1+⋯+b_0c_n=∑_k=0^nb_kc_n-k

Note that 0≤ s_n≤ n. It is useful to visualize this in terms
of the elementary school shifted tabular form:

    [ [      0 c_0b_0 c_0b_1 c_0b_2 c_0b_3      ⋯;               c_1b_0 c_1b_1 c_1b_2      ⋯;                      c_2b_0 c_2b_1      ⋯;      +                      c_3b_0      ⋯ ];                                                                                                                                                                               ;                                                                                                                                                    [ s_0 s_1 s_2 s_3 s_4   ⋯ ] ]

 This makes clear the shift-and-add form. The value of each individual
s_n can be visualized as a stack of blocks. For the special case
of K=0.111⋯=1 one has that s_n+1=∑_k=0^nb_k,
that is, it is simply the total number of one-bits in the first n
locations.

The final step is to reduce the the sum series s_n to a bit-string.
This is accomplished recursively, by performing a carry operation:

    d_n=s_n+⌊d_n+1/2⌋

where ⌊ d⌋ =d1 denotes the floor of
d (the integer part of d). The desired bit sequence is then


    a_n=d_n 2

Equivalently, a_n is the remainder, the part of d_n that
was not propagated to the next location. Explicitly, is is a_n=d_n-2⌊ d_n/2⌋.
The carry-sum propagation can be imagined as a kind of bulldozer,
razing the towers d_n until they are one block high, pushing
the razed bits off to the next location. The resulting sequence (a_n)
is then the bit-string for the product Kx. That is, 

    Kx=0.a_0a_1a_2⋯=∑_n=0^∞a_n2^-n-1


The problem with this algorithm is that the relation <ref>
for the d_n is infinitely recursive, and in general is not guaranteed
to terminate. One has to start at n=∞ and move backwards from
there. There are two plausible scenarios for computing the a_n
in practice. One is to search the n until one finds that spot where
⌊ d_N+1/2⌋ =0; one can then obtain the
a_n for all n<N without issue. The problem here is to find
such an N.

The other way to compute is to observe that the iteration is convergent.
The recursion <ref> only depends on a finite and fixed
number of bits “behind it”, roughly equal to log_2n bits
that come after this. As noted earlier, 0≤ s_n≤ n and likewise,
0≤ d_n≤2n+1. To write down d_n, one needs at most C=1+⌊log_2(2n+1)⌋
bits. This implies that a given d_n can only perturb at most
C-1 bits downstream of it. That is, d_n-C+1 depends on d_n
but d_n-C does not. Thus, in order to correctly compute all bits
a_k for 0≤ k≤ n-C, it is sufficient to set d_n to
some arbitrary value (less than 2n+2) and then iterate (using the
correct values for s_k when k<n). At the end, discard all
d_k and a_k for n-C<k, as they are incorrect.



 §.§ Tree-view


Points:

1) adding one bit is like shifting the tree over sideways.

2) multiplying by one bit is like shifting the tree down-left.

3) adding a number to itself is like shifting tree up (since its just
2x)



§ SIMPLIFIED MODELS OF MULTIPLICATION


The shift-and-add algorithm is obviously rather complex; can it be
replaced by something simpler? The particular question to ask is how
much of the chaotic dynamics of the beta shift is due to the propagation
of the carry bit, and how much of it is due to other parts of the
algorithm? Specifically, the addition of two numbers, which requires
a carry bit, can be replaced by a bit-wise XOR of their bit strings:
this generates “almost” the same results as addition, when the
number of 1-bits in the strings are sparse, but are wrong when 1-bits
appear in the same location: the XOR discards the carry bits. Thus,
a simplified model of multiplication would the the shift-and-XOR model:
it proceeds the same way as shift-and-add, but replaces addition with
XOR. What does this look like, and how does the equivalent of the
beta shift behave under this operation?



 §.§ Shift-and-XOR


The shift-and-XOR algorithm must like the shift-and-add algorithm,
except that it drops the carry bits. Starting from the same spot,
let 0≤ K≤1 and 0≤ x≤1 and represent both by their binary
expansions:

    x=0.b_0b_1b_2⋯=∑_n=0^∞b_n2^-n-1

and


    K=0.c_0c_1c_2⋯=∑_n=0^∞c_n2^-n-1

where the b_n and c_n are either 0 or 1.

Define s_0=0 and s_n+1 to be the result of XOR-ing instead
of adding the bits.

    s_n+1=b_nc_0⊕ b_n-1c_1⊕⋯⊕ b_0c_n=⊕_k=0^nb_kc_n-k

Here, the oplus symbol ⊕ denotes the XOR operation. Note that
each s_n is either zero or one. Reconstructing a real number
from this, one defines

    K⊗ x=0.s_0s_1s_2⋯

where the otimes symbol ⊗ is pressed into service to indicate
the shift-and-XOR product. Note that it is symmetric: K⊗ x=x⊗ K
and so behaves at least a bit like ordinary multiplication. Its is
not distributive over ordinary addition: (a+b)⊗ x a⊗ x+b⊗ x
but it is distributive over XOR: (a⊕ b)⊗ x=(a⊗ x)⊕(b⊗ x).
It is illustrated in figure <ref>.



The range of the shift-and-XOR operation is fundamentally different
from multiplication. First, because the carry bit is dropped, one
has that s_0=0 always, and so that K⊗ x≤1/2 always,
even when both K→1 and x→1. Next, for any value of 1/2<K≤1,
the range of K⊗ x runs over the entire interval [0,1/2]
as x runs over the interval [0,1]. The measure is not compressed
(other than by a factor of 2) , as there is in ordinary multiplication.
That is, if S⊂[0,1] is a measurable subset of
the unit interval, with measure μ(S), then one has
μ(K⊗ S)=μ(S)/2. There are several
ways to prove this. One formal approach is to consider the correspondence
between the natural measure on the reals, and the measure of cylinder
sets on the product topology. That is, the Cantor space { 0,1} ^ω
is endowed with a natural topology, the product topology. The open
sets of this topology are called “cylinder sets”. Their measure
is uniformly distributed over unit interval, precisely because the
Bernoulli shift is ergodic: the one implies the other.

Indeed, the shift-and-XOR algorithm can be best thought of as a formula
for shuffling the bit-strings around, without actually altering them:
re-ordering them, not changing them. The intuitive key to this is
to observe that subtracting x from 1 just re-orders the unit interval,
top to bottom, and that this is the same as flipping all zero bits
to one, and v.v. That is, 1-x=x⊕0.111⋯.

Another way to see this shuffling is to note that a⊕ a=0 and
that 0⊕ x=x. Thus, for a fixed value of a, the string x
and the string a⊕ x are paired together, in a unique way,
so that either can be gotten from the other. The function a⊕[0,1]→[0,1]
sending x↦ a⊕ x is an exchange of these unique pairings
of strings. It is not just a bijection, it is an involution. If the
strings are given their natural lexicographic sort order, the mapping
x↦ a⊕ x is just a certain kind of shuffle of the sort
order; it neither adds new strings, nor deletes any, nor changes their
number. The function is one-to-one and onto. The multiply-and-XOR
algorithm is just a repeated sequence of XOR's: 

    K⊗ x=(c_0x/2)⊕(c_1x/4)⊕(c_2x/8)⊕⋯

and so K⊗ x is nothing more than a reshuffling of strings
(along with a right-shift equal to the number of leading zero-bits
in the binary expansion of K; the right-shift commutes with the
measure on the product topology.) Thus, K⊗ x preserves the
measure on the unit interval (up to a factor of 2^-n due to the
above-mentioned right-shift). That is, for 1/2<K≤1, this discussion
shows that μ(K⊗ S)=μ(S)/2.



 §.§ Self-similarity


There are several self-similarity properties of the shift-XOR worth
noting. It behaves very much like a classic dyadic fractal. Thus,
one has that 


    K⊗(x/2)=1/2(K⊗ x)=1/2K⊗ x

In addition... TODO: illustrate the other symmetry.



 §.§ Similarity Transformations


The shift-and-XOR algorithm acts as a permutation on bit-strings.
As a result, the XOR-analogs of the beta shift and the tent map become
uniformly ergodic, behaving exactly as the Bernoulli shift. The Frobenius-Perron
solution to these is just the uniform distribution, which is featureless.
All of the structure visible in figures <ref>
and <ref> is entirely due to the dynamics of the
carry bit. Effectively, the carry-bit algorithm alters the uniform
distribution of the Bernoulli shift (equivalently, the uniform distribution
associated with the natural measure on Cantor space.) 

Define the XOR-analog of the beta shift as

    c_β(x)=
    2β⊗ x    0≤ x<1/2
    
    2β⊗(x-1/2)    1/2≤ x<1

The factor of 2 makes up for the fact that shift-XOR effectively drops
the top bit; thus the goal is to map each half of the unit interval
into the entire interval [0,1]. 

Given a fixed β, define ⊠_β:[0,1]→[0,1]
as

    ⊠_β(x)=β⊗ x

As observed previously, ⊠_β is an automorphism of
the unit interval, and more: it is a permutation on Cantor space.
Let b(x) be the Bernoulli shift of eqn <ref>;
then one has that c_β=⊠_β∘ b. Taken together,
this implies that the ergodic properties of iterating on c_β
follow directly from the ergodic properties of the Bernoulli shift;
a shuffle, any shuffle on the Cantor set should not alter these ergodic
properties.

TODO: similarity transforms on the transfer operator... and the non-alteration
of the eigenspectrum, even as the eigenfunctions are altered.



 §.§ Multiplication on the Cantor Space


The previous set of results indicates that all of the structure in
the bifurcation diagrams of <ref>
and <ref> is entirely due to the dynamics of the
propagation of the carry sum. To explore this, the notation needs
to be improved on.

The beta shift can be decomposed into multiple distinct stages. First,
there is a conversion from the unit interval to the Cantor space;
this was defined at the very start, but now we need a less awkward
notation for it. Let 

    [        π:2^ω            →        [0,1]; 0.b_0b_1b_2⋯            ↦            x ]

be the projection from the Cantor space to the real-number unit interval,
given by eqn <ref>. Note that it is a surjection:
dyadic rationals (rationals of the form m/2^n) correspond to
two distinct bit strings. For example, 1/2 can be represented as
both 0.1000⋯ and as 0.0111⋯. Cantor space covers the
unit interval. Write the inverse mapping as

    [   π^-1:[0,1]            →          2^ω;            x            ↦ 0.b_0b_1b_2⋯ ]

As a function, it is injective but not surjective. It is usually convenient
to ignore this, and to pretend that both π and π^-1 are
bijections, even though they are not. This rarely leads to practical
difficulties, as long as one stays conceptually tidy. Better yet,
just perform all work on the Cantor space, and project to the unit
interval only when needed.

Next, turn to multiplication. This has three parts. First, the summation
of the carry bits: 

    [         S_β:2^ω               →             ℕ^ω;    0.b_0b_1b_2⋯               ↦ (s_0,s_1,s_2,⋯) ]

where the summation is given by eqn <ref>. Here, ℕ^ω
is Baire space, the space of all infinite-length sequences of non-negative
integers. In number theory, this would be called the space of arithmetic
functions. The second part of multiplication is the propagation of
the carry bits. Denote this as 

    [           C:ℕ^ω               →             ℕ^ω; (s_0,s_1,s_2,⋯)               ↦ (d_0,d_1,d_2,⋯) ]

which is defined in eqn <ref>. Finally, one extracts
the remainder, after propagation:

    [           A:ℕ^ω               →             2^ω; (d_0,d_1,d_2,⋯)               ↦ (a_0,a_1,a_2,⋯) ]

which is given by eqn <ref>. Of the three parts into
which we've decomposed multiplication, only the first part is parameterized
by K. Thus, multiplication, on Cantor space, can be written as
M_β=A∘ C∘ S_β. The shift-and-XOR algorithm
omits the propagation of the carry sum. On Cantor space, it is just
⊠_β=A∘ S_β: the XOR is just modulo-2 of
the carry sum.

To obtain multiplication on the real-number unit interval, we need
merely to re-project from Cantor space to the reals. Thus, multiplication,
given in eqn <ref>, decomposes into

    M_β=π∘ A∘ C∘ S_β∘π^-1

The beta shift of eqn <ref> is then 

    T_β=π∘ A∘ C∘ S_β∘π^-1∘ b

where b is the Bernoulli shift. To simplify notation, it is convenient
to go ahead and provide a symbol for the shift operator:

    [           B:2^ω               →             2^ω; (b_0,b_1,b_2,⋯)               ↦     (b_1,b_2,⋯) ]

so that b=π∘ B∘π^-1. The corresponding beta shift
on the Cantor space is

    B_β=A∘ C∘ S_β∘ B

which eliminates the pesky projection π. It should be clear that
S_β is an injection, the propagation operation C and the
remainder A are both surjections. 

As noted, the shift-and-XOR algorithm can be written as ⊠_β=A∘ S_β;
the step where the carry bits are propagated is dropped. The XOR-version
of the beta shift is


    c_β=⊠_β∘ B=A∘ S_β∘ B

Thus, in this new notation, it reaffirms that B is the true source
of ergodicity, and that A∘ S_β being a permutation does
not alter the basic ergodic property of B. All of the structure
in the bifurcation diagrams can be blamed on the propagation operator
C. 



 §.§ Propagation games


Pinning the “blame” of complex dynamical structure on the propagation
of the carry bits seems to be an open invitation to replace the propagation
operator C by just about anything, to see what happens. Figure
<ref> illustrates some of the things that
can happen. 



Reviewing the images there makes it clear that although fiddling with
the carry bit fundamentally alters point trajectories, it completely
fails to open any doors that would provide insight into the structure
of the transfer operator. The pictures are pretty, but appear to be
meaningless.



§ SCI-FI DAY-DREAMING


This section provides two day-dreams inspired by this material. They
are just that: daydreams. If you don't like fictional daydreaming,
you won't like the material here. Sorry about that.



 §.§ Limits to computation


There are many limits to computation. One limit is the speed of light.
In current generation CPU chips, clock rates in the vicinity of 3
gigahertz=3×10^9 cycles per second. By comparison, the speed
of light in a vacuum is about 3×10^8 meters per second.
Dividing, one finds that light can travel about 3×10^8/3×10^9=10^-1
meters, or about four inches: a bit bigger than the actual physical
dimensions of a chip (typically around half-an-inch on a side), but
not by much. Of course, the speed of light in a metal conductor is
lower – about half the speed in a vacuum. And transistors are small
– more than twenty-thousand times smaller. So, measured in terms of
the size of the transistor, the speed of light is about ten or twenty
transistor-widths per clock-cycle. So, OK, its still fast, at that
length scale. But not really all that fast. The point here is that
the speed of light is a potential limit to the speed of computation,
and it is not all that far away.

In this setting, one can imagine the situation where the speed of
propagating the carry bit during multiplication becomes a limiting
factor. The above work hints at a somewhat boggling idea: can multiplication
be effectively parallelized by working with transfer operators instead?
That is, the multiplication of two numbers corresponds to point-wise
particle dynamics: a discrete particle following a chaotic path through
a complex numerical computation. By contrast, the transfer operator
describes how a distribution propagates through a computation: it
effectively performs “an infinite number” of multiplications at
the same time, in parallel. That is, rather than asking how single
values propagate, one could, and perhaps should, ask how distributions
propagate – parallelize multiplication (for example) to an “infinite”
degree. It is this rather ridiculous idea that suggests that the above
explorations are not purely abstract, but have a potentially practical
application. As I suggested – its a bit of science-fiction day-dreaming
at this point. But it does hint at an alternate model of computation.

Variants of this model have already been explored, for decades. For
example, Crutchfeld defined “geometric state machines” as generalizations
of finite state machines, where, instead of having a finite matrix
(a “transition matrix”) act on a finite vector (the “state vector”),
one instead considers operators acting on homogeneous spaces – that
is, applying a sequence of such operators on homogeneous space. The
most famous and celebrated such space would the ℂℙ^n
– complex projective space, with the operators that act on it being
the the unitary ones: U(n) – such a system defining the n-qubit
quantum state machine. Distributions on ℂℙ^n are mixed
states – and the idea of quantum computing is to evolve such states
through a set of operations. 

The point here is that computation, by means of the time-like evolution
of distributional densities, is already being explored, but in a rather
different context than the one explored here. Here, it seems like
we are bowled over by the complexities of a seemingly much simpler
system.



 §.§ Wave function collapse


There is also a different, bizarrely hypothetical way in which all
of this apparatus could manifest itself. Currently, in order to avoid
the rather severe issues associated with the concept of quantum-mechanical
wave-function collapse, the (vast?) majority of practicing physicists
believe in the many-worlds hypothesis. Clearly, this belief is entirely
correct for microscopic systems, isolated from the usual thermodynamic
hustle and bustle (chlorophyll, rhodopsin and the magnetically sensitive
cryptochromes notwithstanding). But it seems to fly in the face of
daily experience, where we are aware of just one reality. One of my
favorite hypotheses is that this is the result of the (rapid) decay
of macroscopic quantum states down to a probability of zero. The mechanism
is presumably that of decaying subshift measures. Penrose argues that
this has something to do with gravity; but we can go one better: the
natural setting for shift spaces are hyperbolic spaces, as that is
where there is enough room to “fit everything” in a uniform way
consistent with a metric. Curiously, the world we live in – Minkowski
space, is hyperbolic. This suggests that the Many Worlds interpretation
is exactly right, as long as one truly is in Minkowski space, but
that gravitation, which essentially bends or distorts it, squeezes
down the room available for multiple quantum states, effectively forcing
the collapse in this way. 

Put another way: the standard treatment for quantum field theory is
the Feynman functional integral; it can be viewed as an integral over
all possible paths that a “particle” might take. The daydream
is to equate a specific path with the idea of point-dynamics in an
iterated function. As long as one considers only points, and there
movement, one can be completely unaware of either the invariant measure,
or of the decaying eigenstates of the shift operator. In a standard
QFT textbook, all equations appear microscopically time-reversible.
There's almost no idea of a measure, except for the exp-iħ S
in the Feynman integral. The incorporation of gravity into this is
famously difficult. The daydream here is that gravity manifests itself
as eigenfunctions that live off of the shell of unitary evolution.

There is some practical hope of bringing this daydream to fruition:
the theory of subshifts has seen dramatic advances over the last few
decades, getting increasingly abstract, and gaining a firm footing
in very general settings: viz not just in metric spaces, but even
in more general topological vector spaces, and specifically in stereotype
spaces, where most of the devices used in analysis can be exercised
in reasonably safe manner. The point here is that most of QFT can
be formulated using these more-or-less conventional tools and notations.
The trick is to locate and extract those parts that renormalize to
zero, not unlike some of the formally divergent sums explored above,
which can none-the-less be regulated and made to give reasonable answers.
Or at least, that's the daydream. Clearly, got far to got before it
can be reality.



§ TOPOLOGICAL PUSH-FORWARD


The transfer operator is most generally and correctly defined as an
operator acting on the topology of a space, and specifically, as the
push-forward of the (uniform) measure by the iterated function. That
is, given any open set belonging to the topology, the transfer operator
assigns a different open set of the topology: it is a map of sets
to sets. For iterated maps on the unit interval, it is essentially
a map of cylinder sets, the open sets of the product topology. The
shift-XOR experiment shows that the ergodic properties arise from
the Bernoulli shift, and that all other properties, commonly called
“chaotic”, are really the side effect of something else, entirely:
the internal structure of the transfer operator. Fiddling with the
carry-bits cannot reveal this structure; instead, they just define
other, pointlessly goofy iterated functions. Point trajectories fail
to reveal the internal structure of the transfer operator, and at
best point in a misleading direction. To understand the transfer operator,
it must be tackled for what it is: one must look at how intervals
are mapped to intervals, and what sort of symmetries can be discovered
in this mapping. (I've given one sketch of a proof of the transfer
operator as a push-forward in this reference:<cit.>.
There are must surely be better, more accessible and more formal and
mathematically refined presentations; if you, reader, know of such,
please drop me a line.)

The action of the transfer operator on the sets belonging to the topology
of the reals reveals several distinct kinds of actions. The topology
on the reals can be generated from a basis consisting of connected
sets. The transfer operator will map some connected sets to other
connected sets, simply moving them around, shrinking or expanding
them. In other cases, a connected set will be split into two disjoint
parts. For maps that are continuous, there must be regions that have
fixed-points and period-doubling routes to chaos: these correspond
to the (countable number of) “trouble spots” illustrated in section
<ref>.

It seems reasonable to argue that each of these different kinds of
moves creates a distinct group (or monoid) of transformations: in
a certain sense, those transforms that do not change th connectivity,
nor do any folding, are all similar to one-another. It should be possible
to write down exactly which sets belong to this type, and then give
explicit transformation properties between them. Likewise, those connected
sets which are split in two are all similar. It seems like there should
be a prototype: a generic split, followed by some re-arrangement of
the two parts. How can this classification be written in an insightful,
useful way?

I believe that there has been a sufficient number of advances in the
theory of subshifts so that the above vague sketch can be presented
in a fairly concrete way. Unfortunately, most of the relevant material
remains rather arcane and abstract, lacking in direct accessibility
to casual students. I am not currently aware of any adequate yet accessible
treatment.



§ CONCLUSION


What, exactly, is the point of analytic mathematics, especially in
the computational age? Can't one just get a fast computer, iterate
on the logistic map, and find out everything there is to find? Well,
of course, yes, and no: these questions can be taken as either silly
or as deeply philosophical, and it is worth the effort to understand
them and address them properly. 

First, lets dispose of some obvious mis-perceptions. If one carefully
scrutinizes figure <ref>, one
will find signs of a slight unevenness in the horizontal bars. These
are numerical artifacts due to statistical under-sampling: they smooth
out and fade away with additional sampling of the iterated equations.
There is a way to obtain this same figure, far more rapidly, and without
this particular form of numerical noise: one can instead iterate on
equation <ref>. This suggests one philosophical answer:
the goal of mathematics is to find faster ways of computing things;
to discover better algorithms.

A uniting theme between this, and the other text that I have written
on fractal issues, is that they are all explorations of the structure
of the Cantor set, the structure of the space of infinite sequences
of symbols, and the structure of the continuum. That is, we know the
continuum in two different ways: one way is by means of the natural
topology on the real number line; the other is the product topology
on the space of binary strings. The former is suggested by the physical
universe that we actually live in: a continuum with spatial extents.
The latter is suggested by the notion of time and repetition: the
making of choices naturally leads to a tree structure; tree structures
necessarily embed in hyperbolic spaces; the Minkowski space that we
live in is hyperbolic, and this is why, every day, as time passes
on, we get to make new choices precisely because the amount of room
for possibilities is ever-increasing as time flows forward.

What, exactly, do the words “exactly solvable” really mean? So,
for example, equation <ref> involves summation
and multiplication, which has this aura of comfortable preciseness
that an iterated function somehow does not. Where does this sensation
come from? When performing arbitrary-precision numerical computations,
it should be clear that neither addition nor multiplication are simple
or easy: they both require fairly complex algorithms to implement,
and have not-insignificant run-times. To be more precise: the algorithms
are non-trivial because one is using a binary digit expansion to provide
a model for a real number. Different representations of the real numbers
potentially offer different algorithms and performance profiles. One
could represent reals by rationals, but then two issues arise. One
is that the rationals are not evenly distributed across the real number
line: rationals with small denominators cluster about in a fractal
fashion, studied deeply in number theory. As a result, one promptly
gets stuck in a quagmire of trying to understand what a “uniform
distribution” should be. Binary expansions are more “obviously”
uniform. A more basic issue is that, if working with rationals, one
must somehow accomplish the addition or multiplication of two integers.
To accomplish this, one has to represent the integers as sequences
of bits, which only takes us back to where we started. There is no
computational oracle that automatically knows the sum or product of
integers: it has to be computed. The analysis being done in this text
is a kind of a game, where not only is one algebraic arrangement of
symbols is being compared to another, but also one computational algorithm
is being compared for another. Unfortunately, this latter comparison
is very nearly opaque and hidden. If only it could be made visible
in some simple fashion.

The situation here is more easily illustrated in a different domain.
The hyper-geometric series was presented and studied by Gauss; then
Kummer, Pfaff and Euler observed various identities yoking together
different series. By the 1950's, thousands of relations were known,
along with some algorithms that can enumerate infinite series of relations.
The curious situation is that there is no known algorithm that can
enumerate all such relations; there is no systematic way to classify
them. The situation does seem to make clear that there is an interplay
between infinite series and algorithmic relationships between them.
Stated a different way: hyper-geometric series are self-similar, and
the identities relating them are expressions of that self-similarity.

To further sharpen this idea: the dyadic monoid is the generator of
self-symmetry in many common fractals; this is “well-known”, and
I have explored this in other texts. A more general setting for fractal
self-similarities is given by tilings of the hyperbolic surface: to
each tiling, there are corresponding fractals, the self-similarity
of which are given by the tiling. The figures <ref>,
<ref> and <ref>
are clearly self-similar in some obscure way: it is visually clear,
but providing a simple algebraic expression describing the similarity
is difficult; I have not been successful in this. None-the-less, it
seems self-evident that it will be the dyadic monoid that is somehow
responsible for the symmetries, underlying them (unless, of course,
there is some other, as yet undiscovered structure). 

The meta-question is: what is the correct framework by which one can
best understand the interplay between symmetries, infinite series
and algorithms? The current tool-set seems impoverished: it does not
“solve” the systems in this text. Worse, current mathematical
practice reifies addition and multiplication into oracular operations
that magically obtain “the right answer”, when it is clear from
numerical methods that addition and multiplication are necessarily
algorithmic operations performed on finite truncations of infinite
series. It would be nice to place these operations on equal footings,
so as to expose the true nature of this beast.



§ BIBLIOGRAPHY


The references below provide a bibliography that attempts to touch
on all the different ways in which the beta transform and beta expansions
have been studied. Search engines exist to help you find the things
you don't know, and want to find out more about.

tufte


