Here I want to show another argument, why convergence occurs with the type of iteration as in the title.
Only let me rename the involved variables to my own year-long use for my own comfort when writing this:
$b$: I use $b$ for "(b)ase" in $z_{k+1}=b^{z_k}$ with complex numbers $z_k$
$t$: Then I use $t$ for the "fixpoint" such that $t=\lim_{k \to \infty} z_k$
$\qquad $(if this is convergent, otherwise if the inverse iteration $z_{k+1} = \log_b(z_k)$ converges to the fixpoint or at least if the Newton-iteration gives such a fixpoint)
$u$: For the log of the fixpoint I use $u$ such that
$ \qquad b= t^{1/t} = \exp(u \cdot \exp(-u))$ or
$ \qquad t=\exp(-W(-\log(b)))$ and
$ \qquad u=-W(-\log(b))$
conjugacy: Instead of writing $z_{k+1}=b^{z_k}$ it is equivalent to write $y_{k+1}=t^{y_k}-1$ and use conjugacy
$\qquad y_k = z_k/t-1$ and $z_k = (y_k+1)\cdot t$.
$\qquad$ Apart of the advantages which I exploit below, it seems from some heuristics that this might also be numerically more stable than the iterated computation of the un-conjugated original function.
For the conjugated version we can generate the Schroeder-scheme (see for instance wikipedia) for the implementation of the iteration
(which, if $t$ is real and $t \in (e^{-1},e) $ can even be fractionally iterated) . With introduction of the Schroeder-function $\sigma()$ and its inverse $\sigma^{-1}()$ this means to compute for some iteration "height" $h$:
$$ y_{k+h} = \sigma^{-1} (u^h \cdot \sigma (y_k)) \tag 1$$
The original trajectory $\{z_0,z_1,z_2,\cdots \}=\{1,b,b^b, \cdots \}$ (possibly converging to $t$) conjugates to $\{y_0,y_1,y_2,\cdots \}=\{1/t-1,b/t-1, \cdots \}$ (possibly converging to $0$ which is $0 =t/t-1$ by the conjugacy) . Inserting $y_0$ in the Schroeder-equation (1) we have
$$ y_{h} = \sigma^{-1} (u^h \cdot \sigma (y_0) \tag 2$$
The iteration "height" $h$ occurs here only in the exponent of $u$.
The Schroeder-function $\sigma()$ has an invertible formal power-series for $|u| \ne 1$ with constant term $s_0=0$.
Now, if $|u| \lt 1$ then for $h$ going towards positive infinity the cofactor $u^h$ decreases to zero and the whole expression gives
$$ \sigma^{-1} ( 0 ) = 0 \tag 3$$
and conjugacy gives $ (0+1)\cdot t = t$ which is the fixpoint $t$ of the original iteration.
So this should be sufficient for another proof of convergence (according to the OP's title) "on one's own" .
Problems occur, if $|u|=1$.
In the denominators of the coefficients $s_k$ of the Schroeder-function $\sigma(y) = s_1 y/1! + s_2 y^2/2! + s_3 y^3/3! + ... + s_k y^k/k! + ... $ we have products of $(u^i-1)$ for $i=1..k$, so if $u=1$ or some $u^i =1$ all following coefficients become singular and the Schroeder-function $\sigma()$ does not exist for such $u$ where for some $k \in \mathbb N^+$ we have $u^k=1$. (see for more explanation of this my essay functional iteration (pdf) explicitely for the details of the denominators pg 25/26 or more concisely the other essay Eigendecomposition (html) )
However, when $|u|=1$ and $u=\exp(2 \pi î /c)$ with some irrational $c$ the denominators of the coefficients of $\sigma()$ don't show that singularities and further analyses based on the Schroeder-mechanism with the conjugation to $y_k$ might thus be possible this way.
Remark: unfortunately this nonexistence of a $\sigma()$ for $u=\exp(2 \pi î /c)$ with $c$ rational does not allow to analyze by this ansatz the reason, why then also convergence occurs. In case I'll find some operational idea for this I'll add this here, although the OP's question is not at this point.