Convergence properties of $z^{z^{z^{...}}}$ and is it "chaotic" $\DeclareMathOperator{\Arg}{Arg}$
Let $z \in \mathbb{C}.$  Let $b = W(-\ln z)$ where $W$ is the Lambert W Function.  Define the sequence $a_n$ by $a_0 = z$ and $a_{n+1} = {a_0}^{a_n}$ for $n \geq 1$, that is to say $a_n$ is the sequence $z, z^z, z^{z^z} ...$
I am trying to classify the sequence $a_n$ as convergent or divergent (i.e. not convergent) according to $|b|$.  My approach to this involves defining $2$ additional sequences: $b_n = \ln a_n$ and $c_n = b_n - b$.  Then $b_{n+1} =  \ln \left(z^{e^{b_n}}\right) = e^{b_n} \ln z$ and $c_{n+1} = e^{(b + c_n)} \ln z - b= (e^b \ln z) e^{c_n}  - b  = b e^{c_n}  - b = b (e^{c_n}  - 1)$.  Suppose $b_n \to b \implies c_n \to 0 \implies e^{c_n} - 1 \sim c_n \implies c_{n+1} \sim bc_n$.  From here there are $3$ cases to consider:
If $|b| > 1$ then $c_n$ is increasing without bound, so it clearly cannot converge except if $c_0 = 0$ which would imply $c_n = 0$ for every $n$ which is not the case because $b_0 =\ln z$ so $b_1 = e^{\ln z}\ln z \ne b_0$ so $c_1 \ne c_0$, a contradiction $\implies b_n$ is divergent $\implies a_n$ is divergent.
If $|b| < 1$ then $c_n$ is decreasing towards $0 \implies b_n$ is convergent $\implies a_n$ is convergent.
If $|b| = 1$ then $c_n \to c \ne 0 \implies b_n$ is divergent $\implies a_n$ is divergent.
I would like to give credit to user1952009.  The method involving the sequences $b_n$ and $c_n$ is due to them.
Edit 1: After reading Yiannis Galidakis's answer several times, I have decided it is mostly satisfactory.  Having said that, I feel I must clarify my position on chaotic sequences.  When I said the concept of a chaotic sequence was nonsense, that was an exaggeration.  What I meant to say was the definition given in the previous post didn't really make sense to me.  The one given here is detailed and precise enough that I can understand what it actually means.
However, I do not believe the statements "$a_n$ is chaotic" and "$a_n$ does not converge" are incompatible.  Setting aside the question of whether $a_n$ is chaotic or not, I feel I have given a plausibility argument (although not a rigorous proof) that $a_n$ is convergent if and only if $c_n \to 0$ which happens if and only if $|b| < 1$.
On a different note, I have been experimenting numerically with $a_n$ for at least $3$ years.  I believe (but cannot prove) that whenever $z$ is not purely real and $|b| = |W(-\ln z)| > 1$ the sequence $a_n$ 
actually has a $k$-cycle for some $k \in \mathbb{N}$.  It appears that $k = 3$ for all sufficiently large $|z|$ and that $k \to \infty$ as $|b| \to 1$.
Similarly, whenever the sequence "shoots off" too a neighborhood of $\infty$ I have found that $a_n \sim$ {$..., \infty, 0, 1, z, z^z ...$}  I will admit however that I cannot predict when the sequence will "shoot of" to $\infty$.
At this point I have $2$ main questions: is my convergence analysis correct?  Are there any values of $z$, not purely real, such that $a_n$ is unbounded, i.e, the sequence "blows up" to $\infty$ and never "comes back down."
Apologies for a somewhat lengthy post.
Edit 2: In Addendum 2 in Yiannis Galidakis's answer, it is stated that if $c \in \mathbb{C}$ and $|\Re(c)| > \exp(\exp(-1))$ and $\Arg{c}=2\pi\alpha$ with $\alpha\in\mathbb{R}\setminus\mathbb{Q}$ then the iterated exponential will be unbounded.  This is not quite correct.  The part about irrational rotation is fine, but it is not enough to have $|\Re(c)| > \exp(\exp(-1))$.  Consider $c = 2\exp(e\pi i/14) \approx 1.6393177+1.1457037i$.  Clearly $c$ satisfies both conditions. However $|W(-\ln c)| \approx .890512$.  Hence the iterated exponential is convergent, hence it is bounded.  A few other counterexamples are given by $c = k\exp(e\pi i/14)$ for $k \in \mathbb{R}$ and $1.76253 \le k \le 2.34896$
 A: $\DeclareMathOperator{\Arg}{Arg}$
Let me try to explain what happens with these sequences using a simpler example. Modulo the magnitude of the complex numbers they shuffle, they are almost identical behaviorally to the sequence $a_n=\exp(i n)$, $n\in \mathbb{N}$, for example. Let's see that sequence on the complex plane for say, $N=250$.
restart;
with(plots);
a := proc (n) options operator, arrow; exp(I*n) end proc
N:=250;
L1 := [seq([evalf(Re(a(n))), evalf(Im(a(n)))], n = 1 .. N)]
plot(L1, style = point, symbol = cross, scaling = constrained)

This sequence is dense on the unit circle. It is an iterated rotation with angles which are irrational multiples of whatever unit angle you use. Accordingly, their projections, $\Re(a_n)$ and $\Im(a_n)$ are dense in their corresponding intervals, $[-1,1]$ and $[-i,i]$. These are standard results which you can find in this forum, on the question "does $\sin(n)$ converge", etc.

So the above sequence, is afortiori divergent, using your favorite catchall for any non-convergent sequence.
This behavior is not the problem though with the iterated exponential. The difference between the behavior of $a_n$, above and those other sequences you see on the iteration of the exponential map, is more specific.
Suppose then that we generate additionally a unit impulse along the direction of the angle $\theta=\pi/3$. Here's the corresponding impulse:
r := proc (x, epsilon) options operator, arrow;
piecewise((1/3)*Pi-(1/2)*epsilon < x
and x < (1/3)*Pi+(1/2)*epsilon, 1+1/epsilon, 1) end proc
plot(r(x, .9), x = -Pi .. Pi)


Now, let's see the corresponding distribution for the sequence $a_n$:
L2 := [seq([evalf(Re(r(.9, n)*a(n))), evalf(Im(r(.9, n)*a(n)))],
n = 1 .. N)];
plot(L2, style = point, symbol = cross, scaling = constrained);


You can easily verify now, that the sequence has acquired some local "fuzz", with several elements sling-shot to farther values. Now, if you turn the fuzz into locations, by using the option of line-connect on the graph, you will see the overall momentum of the distribution, which is something like this:
plot(L2, style = line, scaling = constrained)


It is fairly clear now, that the momentum "tends" to blow the orbits away. The total trajectory then, can be used to find bounds for the domain of the sequence. In fact, the Fatou domain of the iteration, is determined precisely by the Convex-hull of the above figure (which is symmetric with respect to the axis $0$-farthest point.
Now, if you make your impulse at $\theta=\pi/3$ stronger, say, by picking a smaller $\epsilon$ (because the argument is dense around $\theta=\pi/3$), you can always find points outside the old hull, with a tendency to be moving outwards in the direction of $\theta=\pi/3$.
Obviously then, if your impulse is very strong, you are guaranteed finding some $a_n$, close to some neighborhood of infinity.
Now, make the parallel with the iteration of the exponential and replace the above distribution by the actual distribution, which is none other than the magnitude of the previous term (or some function of it thereof (*)).
The whole story then, is precisely that: When you iterate the exponential of an irrational base, you are basically rotating at irrational angles and at the same time impulsing as strongly as the magnitude of your previous term in some direction $\theta_0$.
Therefore there's a probability $\epsilon>0$ you may be sling-shot to a neighborhood of infinity. But even up to here, the sequence cannot be characterized as "chaotic".
It can be labelled "chaotic", if you may allow the term, when you have wittnessed at least one sling-shot in the neighborhood of infinity. Because when that happens, you cannot tell what your next term will be! In fact, the iteration of the exponential with an irrational winding number, may after that shot, return the sequence to its previous bounded Fatou domain.
The characterization "chaotic" therefore is qualified only as your inability to predict when your sequence will jump close to a neighborhood of infinity and come back.
On one hand you have a sequence which is bounded by $M_k$ in some Fatou domain (which is appropriately called a "wandering Fatou domain" because the sequence neither converges nor is unbounded there) and on the other it may suddently jump close to complex infinity and stay there or return after one iteration.
In particular, it becomes very difficult to predict the sequence:
$$B_n=\{|a_n|,n\in\mathbb{N}\}=\{M_1,M_2,\infty,M_3,\infty,\infty,\cdots\}\text{    (****)}$$
That's what is usually meant as "chaotic" in the dynamics of the iterated exponential, in the sense of that changing the $\epsilon$ of the underlying distribution, may result in humongous changes to the magnitude of $a_n$ (sling-shot the sequence to infinity and back).
Hence the characterization of these domains as "wandering" in the dynamics of the iterated exponential.

Notes:
You can have even worse cases. Under certain circumstances, the Convex-hull of the domain may explode (via a Knaster explosion) and reach all the way to complex infinity consisting of strange indecomposable continua, which are entire regions of positive measure filled with points from the sequence (**).
So, if you still object to the characterization "chaotic" as nonsense, as you say, simply call them "divergent" and let's agree on that (***). Ping me if any of the above is not clear and I will stuff in further explanations.
(*) Using the notation of the iterated exponential by Devaney, as $E_{\lambda}(z)=\lambda\cdot \exp(z)$, for example, this is nothing other than the quantity $\lambda=\ln(c)$, the logarirthm of the base of the iteration.
(**) This happens easily in any neighborhood of the point $c=\exp(\exp(-1))$, for example, with/where $|\lambda|>1/e$.
(***) The gist of the problem being really the inability to apply the theorem: Every bounded sequence contains a convergent subsequence, to extract a convergent subsequence. So you cannot pin-point the nature of such sequences, since you already have problems determining whether they are (eventuallly) bounded or not, by (sequence (****) in text).

Addendum (in context)
Now you can put the above result into context with your original question. The multiplier of the iterated exponential is given by:
$$t=-W(-\ln(c))$$
As you well found out then, this splits into three cases:


*

*$|t|>1$, in which case the iterated exponential diverges (by fixed point iteration)

*$|t|<1$, in which case it converges (to $\frac{W(-\ln(c))}{-\ln(c)}$ (Corless) again by fixed point iteration)

*$|t|=1$, (boundary of Shell-Thron region) which splits into the two cases below:
3a. $|t|=1$ and $t^n=1$, i.e. rotating by rational multiples of $2\pi$, which is taken care by Baker and Rippon (i.e. converges as in 2)
3b. $|t|=1$ and $t^n\neq 1$, i.e. irrational rotation, which gives the analysis above.
Now that all possible cases have been classified, suppose that your advisor asks you the following question:
So, what does the iterated exponential do in case 3b? Pick the right answer from the items below and we are done.


*

*It diverges

*We don't know

*It's chaotic (maybe nonsense by your standards)


Give it a shot. Any answer is fine by me.

Addendum 2 (for your new question)

At this point I have 2 main questions: is my convergence analysis
  correct?

Insofar that it reduces to the 3 cases shown above, yes, it is.

Are there any values of $z$, not purely real, such that $a_n$ is
  unbounded, i.e, the sequence "blows up" to $\infty$ and never "comes
  back down."

Addendum 2a (Addressing your concern in your latest comment)
WhoOps! In response to your question, above, I gave the following answer, which is, nicely wrong:

Trivially. Just pick any $c$ for which $|\Re(c)|\gt \exp(\exp(-1))$
  and $\Arg(c)=2\pi\alpha$, with
  $\alpha\in\mathbb{R}\setminus\mathbb{Q}$. For such a $c$, the iterated
  exponential will generate an iterated irrational rotation under an
  unbounded distribution (because $|\Re(c)|\gt \exp(\exp(-1))$, which is
  the upper bound of convergence by Euler), which will eventually
  sling-shot everything to $\infty$ in a spiral manner, by the very
  analysis that I gave you above.
If you saw it on top of the Riemann sphere, it would be a spiral
  converging to the top point of the sphere.}

I answered too hastingly and without checking it first. The condition is obviously wrong, since the strip $e^{-e}\le\Re(c)\le e^{1/e}$ intersects the Shell-Thron region (figure below).

Besides, it's the very nature of those chaotic sequences that you cannot predict what they do. Let me change my answer to the following:
Unknown. As a "general" rule, I'd propose the following: Let $\phi(z)=\exp(z/\exp(z))$, the potential map of the iterated exponential. Then if $t(\rho,\theta)=\rho\cdot e^{i\theta\pi}$, with $\rho>1$, $\theta\in\mathbb{R}\setminus\mathbb{Q}$ and $|\lambda|=|\ln(c)|\le 1/e$ the iterated exponential with base  $c=\phi(t)$ may contain a subsequence which shoots to infinity.
That's as close as I can come with an answer and obviously you have to move away from the unit circle a bit, with $\rho>1$ to get some visible results. I am also excluding a neighborhood $B(e^{1/e},\delta)$, $\delta>0$ (where it may happen $|\lambda|=|\ln(c)|> 1/e$), to avoid dealing with the point $\eta=e^{1/e}$, around which all sorts of weird things can happen (because $\phi$ is not conformal there, so we have Devaney's Knaster explosions, etc).
That's then only a partial answer to your question (for values of $c$ whose iterated exponential shoots to infinity at least once) and not an answer for when the entire thing is blown to infinity.
For the whole thing to go to infinity, I would look at near the far $x$-axis, say at values like $c=M+\epsilon\cdot i$, for large $M$ and various $\epsilon$'s, at both rational and irrational angles, but calculations don't allow me to see anything specific, such as a sufficient condition.
Perhaps you can lend a hand here and investigate it yourself. It would make for a nice result.

Addendum #3: (In response to your answer)
The exact Theorem of Baker and Rippon reads:

$\{{^n}c\}$, $n\in\mathbb{N}$, converges for
  $\lambda=\ln(c)\in\{t\exp(-t)\colon |t|<1\text{, or }t^n=1 \text{, for some }
 n\in \mathbb{N}\}$, and it diverges elsewhere.

It's a result dated 1986, I believe. From the above quote, it is obvious that the authors characterize case 3b) as "divergent", as well. Which isn't of course wrong...but,
since then, many (potentially) chaotic maps (including the exponential of base $\lambda=\ln(c)$) have been analyzed in detail by the latest dad of Complex Dynamics, Robert Devaney, the most prolific author in this field, who also gives exact definitions for when such states may occur.
Hence, the dichotomy between "divergent" and "chaotic" in this case (3b) is carried over only as an artifact of the fact that Devaney's work is more recent to Baker and Rippon's result.
Your frustration and unwillingness to label something as "chaotic" is sort of understandable, in a way. It is a term many people abuse without knowing exactly what's going on, but in this case it is fairly well documented with definitions, etc., from Devaney's work.
Besides, such points exist in much simpler maps, such as in the Mandelbrot map, etc.
A: This is not a new answer, but only intended to give some illustration for the cases, that $|b|=1$ and the two subcases, that 1) $b$ is a rational-order complex root of the unit $b = \exp( 2 \pi î /q) $ where $q \in \mathbb Q$ and 2) $b$ is an irrational-order complex root of the unit.               
Remark: the examples are computed using Pari/GP with internal precision of 200 decial digits and for verification re-computed with 400 digits precision. 
Update: I've extended this discussion 9'2018 and explained it more in this 8-page essay: equator .

Case 1) As Yiannis Galidakis has already pointed out, in the case 1) (where $q \in \mathbb Q$) we'll have convergence (different from the assumption in the OP's question - and also from an earlier answer of myself, in which I suggested non-convergence); two numerical examples suggest that by visual inspection and where the numerical errors also seem to be neglectable.         

The first picture shows the orbit for $q=5$ beginning at $a_0=z=1$ giving a five-edge-star contracting along with the iterations. To see the effect of the contraction clearer I connected each fifth iteration with lines and get this way five separated curves (call them "partial orbit"), each one approximating a/the fixpoint $t_5$ in the middle which can be found numerically using Newton-iteration or using the Lambert-W-formula and setting $t=\exp(u)=\exp(-b)$              
 
However, the rate of contraction seems to diminuish and by visual inspection only this might possibly get stuck in 5 periodic accumulation points.      
Thus I've prepared another picture, where the numbers of iterations per documentation in the curves increase with the powers of 2, so the first round of 5 iterations begin at $a_0$, then the next round at $a_10$, next $a_{20}$, $a_{40}$,$a_{80}$,...,$a_{10240}$... Now we see the contracting effect again, but still with some diminuishing rate:


It is similar with $q=3$ (so $b$ is the 3'rd complex unit-root) giving an orbit with three partial curves, apparently contracting to the fixpoint $t_3$. As before, the number of iterations documented by dots was increased by increasing powers of 2 to focus the general tendency:       
 

Case 2) While the two examples with rational $q$ give clear images with $q$ partial curves, in the case 2) selection of irrational $q$ give more irregular orbits without (so far) discernable contraction or expansion. Instead, when iterating the orbit seems to fill out some curve with fractal border but with a possibly meaningful smooth interpolation in most regions. Periodicity means here, that with a numbers $r$ of steps the points $a_n$ and $a_{n+r}$ are near, and the nearness can be made arbitrarily small by using $r$ from appropriate convergents from the continued fraction of $q$. So with $q=\exp(1)$ I found a nice value for $r=87$ giving the following image, consisting of $r$ partial orbits:
(Remark: the following pictures were computed by shifting the fixpoint to the origin, which also means to use $d_{n+1} = t^{d_n}-1$ instead of $a_{n+1}=z^{a_n}$. The shape of the curve is only marginally affected)
 
The different colors of the lines mark different partial orbits (only a couple of them are shown). I did not discover contraction or expansion of the interpolation-curve so far; also in the top-right region there might be a singularity - but I've not done enough iterations to see this clearer.

An artistic looking image give the value $q=\pi$. Here I did separation of the orbit into 11 partial orbits, and to declutter the image I only show one partial orbit connected by a line.       

Separation into 355 partial orbits, btw., gives very dense partial orbits, but I did not yet produce a visual more satisfying image with this option.
A: This is an extended response to Yiannis Galidakis's question.
I have been experimenting with the value $c$ they gave in the previous post.  It is in case 3b since $|t| = |W(-\ln c)| = 1$ and there is no $n \in \mathbb{N}$ such that $t^n = 1$.  I have evaluated the sequence $a_n$ for $0 \leq n \leq 10^8$.  As in Gottfried Helm's analysis, I have found no evidence that the terms are contracting towards some limit.  In fact, it appears the sequence is frequently "slingshot" to a neighborhood of $\infty$, and that the terms are being "slingshot" progressively farther out as $n \to \infty$.  As an example, the largest such "slingshot" has $|a_n| = 17.1$ for $n \leq 10^6$.  For $n \leq 10^8$  this increases to about $25$.  In all such cases, I have observed that $a_n \sim$ {$...,\infty, 0, 1, z, z^z, ...$}.  For these reasons I believe $a_n$ is divergent (i.e. not convergent) in case 3b.
As to whether $a_n$ is chaotic or not, I am not completely sure.  This may be a bad analogy, but, the way I have come to think of it, when we say a sequence or other mathematical object/phenomenon is chaotic, what we're really saying is the condition number of some function is very high relative to some other function.  By this definition, I believe the sequence $a_n$ can be chaotic in all cases, but is especially chaotic in cases 3a and 3b.  For one thing, an arbitrarily small perturbation in $a_0 = z$ can leave the sequence in the same case or shift it into any of the others.  Also, using the formula $cond(f(x)) =  \frac{xf'(x)}{f(x)}$ we see that $a_n$ is very ill-conditioned when it is "slingshot" near $\infty$.
In spite of these issues, I do not lose much sleep over it.  As far as I know, there are $4$ ways to alleviate the difficulty of an ill-conditioned problem.


*

*Supercomputers

*Arbitrary precision calculators

*A very sophisticated algorithm

*A lot of time and effort


Of course there are some problems where even all $4$ of these are not enough.  As someone once said "Computers and calculators are not replacements for mathematical thought. They are just replacements for some kinds of mathematical labor, either numerical or symbolic. There are, and always will be, mathematical problems that can’t be solved by a calculator or computer, regardless of its size and speed."
