(The presentation that follows is a very light introduction. For more austere results, consult Paper 1).

Consider the relation:

(1)

Replace now the exponent, (2) with the left side of (1). We get:

By repeating the process of this substitution of the last exponent recursively, we obtain finite sequences of n radicals, like:

which are valid for all n > 0.

Let us see what happens "in the limit".

Claim:

In fact, we can prove a much stronger result:

[x^(1/x)]^{[x(1/x)][x(1/x)]...} =

{x, if x is in [1/e,e],
y, where y is in [1/e, e] and y satisfies: y^(1/y)=x^(1/x), if x is in (e,+∞)} (2).

The correct way to "see" this, is to grab the recurrence that describes the above substitution.

Let us define:

a_n=[z^(1/z)]^a_n-1, with initial condition: a₁=z, z in C, z<>0, (3)

We have to be careful with what we mean by (3). To clarify the operations performed, we agree that (3) means:

a_n=e^{{Log(z)/z*a_n-1}}, with initial condition: a₁=z, z in C, z<>0, (4)

where it is understood that we are working with the principal branch (Log(z)) of log(z), on every recursive step.

In view of (4), it is still true that lim_n->+∞a_n=z, as the sequence is again constant and always equal to z.

Since in our recursive definition of (3), z was arbitrary (but non-zero), we can define a_n using z=z_o, for any specific z_o<>0. In particular, a_n is a constant sequence (and equal to z), if we pick a z that satisfies:
z^(1/z)=c, c in C, c<>0, (5)

where it is understood that what we mean by (5) is equivalent to:
e^(Log(z)/z)=c,

and that we are working again with the principal part of log(z), Log(z).

Equation (5) has infinitely many solutions, given by:

-W(k,-Log(c))/Log(c), (6)

where k is an integer, denoting the Lambert W branch, and where it is understood that we are working with the principal branch of log, Log, and W(0,z)=W(z) is the principal branch of the W function.

Here is what Lambert's W function looks like: The principal branch and the pair of branches W(-1,x) and W(1,x) share an order 2 branch point at -1/e. The branch cut dividing these branches is the subset of the real line from -infinity to -1/e, and the values of the branches of W on this branch cut are assigned using the rule of counter-clockwise continuity around the branch point. This means that W(x) is real-valued for x in the range -1/e .. infinity, while the image of -infinity .. -1/e under W(x) is the curve -y cot(y) + y*i, for y in 0 .. Pi.

We now define a new symbol, called the infinite exponential c^ccc..., for c in C, c<>0, as follows:

c*= c^ccc...=lim_n->+∞a_n with a_n as in (3) and with z being the principal solution z_o, of (5): z^(1/z)=c. (7)

It can now be verified that the recurrence (3) with z=z_o, will be a constant sequence, always equal to:

z_o=-W(-Log(c))/Log(c).

Therefore, its limit will clearly be z_o.

The definition is unambiguous, provided we are working with the principal branch of log, Log(z) and with the principal branch of W. For, if c₁=c₂, c₁ and c₂ in C, c₁, c₂ <>0, then -W(-Log(c₁))/Log(c₁)=-W(-Log(c₂))/Log(c₂) which forces c₁*=c₂*.

What we really have done here, is to define the infinite exponential as a multiple-valued function f:C->C, with f satisfying:

f(c)=-W(k,-Log(c))/Log(c), k in Z. (8)

This definition is consistent with the "usual" algebraic infinite exponential, WHEN this exponential converges. For example, the "usual" infinite algebraic exponential 5^55... is meaningless (however we assign it a complex value under the new definition), as it diverges, where as the exponentials (2) for example, converge. (Verify for example, that in the case of (z_o=2^(1/2)), -W(-Log(2^(1/2)))/Log(2^(1/2))=2)

f above in (8), can be then alternatively written as: e^-W(-Log(c)), based on the definition of Lambert's W function:

W(z)*e^W(z)=z. (9)
for z=-Log(c).

(Details about this function (which is subsequently renamed to F(x)) are presented in The analytic continuation of the Hyper Tower function F(x) = ^+∞x).

The function is real valued for real values of z in the interval (0,e^(1/e)], {in fact for [0,e^(1/e)], since lim_z->0f(z)=0}, which gives a good hint for which algebraic infinite exponentials converge. In particular, an algebraic infinite exponential converges, iff its base belongs to the interval [(1/e)^e,e^(1/e)]=[0.06598803588,1.444667861]. (That's why the infinite exponentials (2) above with base x^(1/x) converge. See Lemma below). The largest such real convergent exponential (in absolute value), will thus be the one with base e^(1/e), so it follows that [e^(1/e)]^{[e(1/e)][e(1/e)]...}=f(e^(1/e))=e. (Note that this is also the largest real base an infinite convergent exponential can have, since the maximum of the function x^(1/x) is attained at x=e.)

The smallest such real convergent exponential will be the one with base e^-e.
And in particular [(1/e)^e]^{[(1/e)e][(1/e)e]...}=f((1/e)^e)=1/e.  To see why these are the endpoints of the convergence, we consider the function whose iterates we investigate, i.e. g(z)=c^z. Then, we know that if c_o is a fixed point, then if |g'(c_o)|>1 the fixed point is repellant, and if |g'(c_o)|<1, the fixed point is attractive. Since our fixed points are given by (6), all we need to do is solve the equation |x^{{-W(-log(x))/log(x)}}*log(x)|=1 to find the points for which convergence is guaranteed. Solving this gives x={(1/e)^e,e^(1/e)}.

What happens on the interval (0,(1/e)^e)?

Let's look at some preliminaries:
Consider a small number d close to zero. It is immediately clear that if d is close to zero, d^d will be close to 1. It is also clear that d^dd will be again close to zero, as d^d is close to 1. It becomes apparent thus, that there is some major alternation of the terms between two values close to 0 and 1 of the sequence d^dddd... Can we say anything about those numbers? Yes we can. Consider the two sequences, of EVEN and ODD terms of d^d...

The first sequence will be the sequence an of even number of terms in the exponential:
a_n={d^d, d^ddd,...}
The second sequence will be the sequence of odd number of terms in the exponential:
b_n={d, d^dd,...}

It follows that *IF* a_n converges, its limit will be solution to the equation:
lim1=d^dlim1 (10)

and accordingly, *IF* bn converges, its limit will be the solution to the equation:
lim2=d^lim1 (11)

Now *IF* the infinite exponential is to converge, those two limits must be equal, so:

lim1=lim2 (12)
This means that d^dlim=d^lim (13)

Solving (13) for lim we get

lim=-W(-log(d))/log(d) (14)

This shows that *IF* the infinite exponential converges, its limit must indeed be given by (6) (or (14)). If the exponential does not converge, equations (10) and (11) give the two separate limits of the two cycle.

(Should you wish to see how the equation c^cx=x is solved, please consult the article: Solving the Second Auxiliary Real Exponential Equation).

(Details about the actual Hopf bifurcation that occurs in (0,e^-e) are presented in The analytic continuation of the Hyper Tower function F(x) = ^+∞x).

Lemma:
f(w^(1/w))=w, if w is a real number in [1/e,e].
Proof:
Let z_o be the solution to the REAL equation z^(1/z)=w^(1/w) with z, w real in [1/e,e]. As the function y=x^(1/x) is 1-1 and onto the range of [1/e,e] (see below)=>z_o=w. But z_o=f(w^(1/w)), from which it follows that f(w^(1/w))=w.

Now we look at what happens when the base is of the form x^(1/x) and x ranges in (e,+∞). For this, we have to take a look at the graph of y=x^(1/x). Since y attains a global max at x=e and since y is continuous in the range (1,+∞), and since y=1 is an asymptote as x->+∞, it follows that if x>e, there always exists a w in (1,e], such that w^(1/w)=x^(1/x). THIS w is precisely the limit of our infinite exponential. (An example is given below for 6.62^(1/6.62) which equals 1.330444232=1.562011537^{(1/1.562011537)}.) Verify then, that f(w^(1/w))=f(x^(1/x)), but according to the Lemma above, f(w^(1/w))=w, since w is in (1,e]. And the limit is clearly shown.

In particular, many infinite exponentials can be investigated this way. What is i^iiii... (if anything)? The answer can be determined if we could somehow use the analysis above to describe the exponential. And we can do that. All we have to do is solve i=z^(1/z), which gives z_o=2*i*W(-i*Pi/2)/Pi. In particular then, f(i)= z=z_o=2*i*W(-i*Pi/2)/Pi. (thanks to John Bailey for verifying the solutions). So, i* is essentially a multiply valued function, and its principle value is given by:

i^ii...=2*i*W(-i*Pi/2)/Pi=.4382829366+.3605924718*i. (In Maple, the rest of the "solutions" can be gotten as 2*i*LambertW(k,-i*Pi/2)/Pi, k=-∞..+∞, for the rest of the branches of Lambert's W function).

Following the above discussion, it is natural to ask if infinite exponentiation can produce an interesting fractal. Indeed it does produce beauty, since complex numbers are involved. Upon the suggestion of John Bailey, I set out to investigate possible algorithms which can show the speed of convergence.

My original code was a short Maple excerpt, and it did produce the following image:

The code that produced that was:
>with (plots):
>F:=proc(x,y)
>local c,z,m,f,l,eps;
>z:=0;m:=0;eps:=1e-4;
>c:=evalf(x+y*I);
>if abs(c)<eps then
>m;
>else
>l:=log(c);
>f:=evalf(-LambertW(-l)/l);
> to 100 while abs(z-f)>eps do
>z:=evalf(exp(l*z));
> m:=m+1;
>od;
> m;
>fi;
>end:
> plot3d(F, -1 .. 1, -1 .. 1, grid=[100,100], shading=ZHUE, style=PATCHNOGRID, orientation=[0,0], axes=BOX);

Note that the fixed points for the orbit of c are given by  (6). The above code simply checks how close the iterate of the function f(z)=c^z is to the fixed point.

However the above code is very slow so it is not recommended you run it on Maple.

Robert Israel optimized the above code, by avoiding the Lambert's W function, and by eliminating complexes altogether.
His Maple code is longer but much faster overall.
Here it is:

>with (plots):iterf:=proc(lre,lim)
> local eps,theta,R,m,zre,zim,fre,fim;
> eps:=1e-4;m:=0;
> zre:=0;zim:=0;fre:=1;fim:=0;
> to 100 while (abs(zre-fre)+abs(zim-fim)>eps) and (abs(zre)+abs(zim)<1e37) do
> fre:= zre; fim:= zim;R:=exp(lre*zre-lim*zim);
> theta:=lre*zim+lim*zre;
> zre:=R*cos(theta);
> zim:=R*sin(theta);
> m:=m+1;
> od;
> if abs(zre)+abs(zim)<1e37 then 120-m
> else m-120 fi;
> end:
> F:=proc(x,y)
> local c,z,l,lre,lim,eps,zr,zi;
> z:=0;eps:=1e-4;
> c:=evalf(x+y*I);
> if abs(c)<eps then 0
> else
> l:=log(c);
> lre:=Re(l);lim:=Im(l);
> evalhf(iterf(lre,lim));
> fi
> end:
> plot3d(F,-1..1,-1..1,grid=[100,100],shading=ZHUE,style=PATCHNOGRID,orientation=[0,0],axes=BOX);

The above code produced the same fractal but much faster. (Settings are -4..4,-4..4 and grid=[300,300].) The dark pixels in all fractals below are actually unknown areas.

Finally, I translated Robert's code to Pascal, which produced the fastest results and is a bit more accurate than the Maple code.

Note that the central bulb's real range is [(1/e)^e,e^(1/e)]. The actual reason for this is explained in the next article. This section is shown in detail below.

To see full blown pictures of the above "tetration" fractal generated using FractInt, visit Daniel Geisler's website.

For more mathematical information on this fractal, visit this page.

The above process, is a "Mandelbrot" process.We can instead alter the process to produce Julia sets as well. (i.e. fix c and iterate z taken from the grid, through the same function. The following picture, is the Julia set for c=i. The attractor is located at .4382829366+.3605924718*i. The principal value of i^ii... = -W(-log(i))/log(i). (The little white dot marks the origin).

Here is a better looking image for the above Julia set, produced with SpangFract, a fractal program for the Macintosh, by David Spangler.

Back to Mathematics