In a previous mathstack question I posted a formal series solution for finding the fixed points of $z \mapsto \exp(z)-1+k$. This analytic Taylor series in $\sqrt{-2k}$ can generate both fixed points for bases between $1 < a < \exp(1/e)$ as well as the complex conjugate pair of fixed points for bases $>\exp(1/e)$, and works for complex bases, as well as real bases<1. Since that post, I I have a better way to formally calculate that series.
For the Op's question, it seemed there might be a series that could work in the neighborhood of $\exp(-e)$, to generate the two periodic solutions. I was able to generate just such a series which work wells for $0<a\approx<0.6$. For calculating Gottfried's fixed point pair solution for $a^{a^x}$ for $0<a<e^{-e}$. we need this new Taylor series, which is centered on the parabolic fixed point with multiplier -1. We expect that the solution would be analytic in $\sqrt{x-\exp(-e)}$. But just as with the earlier problem, it is simplest to find such a Taylor series using a mathematically "congruent" problem.
For the Ops problem, $a^{a^x}$, for $0<a<1$ I use this congruence equation, where we want period2 fixed points of $f(y)$ with $$k=\ln(-\ln(a))-1$$ and then we can instead iterate:
$$y \mapsto f(y,k);\;\; \text{where} \;\;\;f(y,k)=-\exp(y)+1+k;\;\;\;\;\text{and}\;\;\;\;y=z\ln(a)+\ln(-\ln(a));$$
$$z \mapsto a^z;\;\;\;\text{is congruent to}\;\;\; y \mapsto f(y,k);\;\;\;\; k=\ln(-\ln(a))-1; $$
$$ f\big(z\ln(a)+\ln(-\ln(a))\big)=(a^z)\ln(a)+\ln(-\ln(a));\;\;\; f \; \text{is congruent to}\; a^z$$
k=0 corresponds to $a=\exp(-e)$, which has a multiplier of -1. The two periodic fixed points of $f(y,k)$ has a series which I shall call $g$, whose definition is below. We then find a formal series solution for the two periodic fixed point pairs of $f(f(x,k))=x$
$$f(x,k) = -\exp(x)+1+k$$
$$f(g(\sqrt{6k},k) = g(-\sqrt{6k})\;\;\;\;\text{g(-x) is the other fixed point}$$
$$-\exp((g(\sqrt{6k}))+1+k = g(-\sqrt{6k})\;\;\;\;\text{definition of f}$$
$$-\exp(g(x))+1+\frac{x^2}{6}=g(-x)\;\;\;\;\text{by substituting k=x^2/6}$$
That is the formal series definition for g, where the $x^2/6$ term was chosen so that the for the g Taylor series below, the x^1 coefficient=1. In this equation, $g(x)$ corresponds to one fixed point, and $g(-x)$ corresponds to the other fixed point.
And the two cycle fixed point of $y\mapsto -\exp(y)+1+k$ may be found by $y=g(\pm\sqrt{6k})$. And then the fixed point pair of z for
$z\mapsto a^{a^z}$ may be found by using this equation:
$$z=\frac{g(\pm\sqrt{6k})-\ln(-\ln(a))}{\ln(a)};\;\;\;\; k=\ln(-\ln(a))+1 $$
The first 16 terms of the series for g are as follows. I wrote a pari=gp program to calculate the formal series for g which requires iterating solving 2x2 simultaneous equations for pairs of consecutive terms, but it is not too bad. With enough terms, the series can be used on its own, or the series can be used as an input to Newton's method to get a more accurate answer.
g=
+x^ 1* 1
+x^ 2* -1/6
+x^ 3* 1/20
+x^ 4* -1/90
+x^ 5* 523/151200
+x^ 6* -23/28350
+x^ 7* 239/1008000
+x^ 8* -19/340200
+x^ 9* 1471949/100590336000
+x^10* -6583/1964655000
+x^11* 94891697/130767436800000
+x^12* -49909/328378050000
+x^13* 18670028801/988601822208000000
+x^14* -520019/241357866750000
+x^15* -88448773393/67224923910144000000
+x^16* 254033333/492370048170000000 ....
So lets say we want to two real fixed points for $a^{a^z}$ for a=0.04
Then $k=\ln(-\ln(0.04))-1\approx 0.1690321758870$, and we want $g(\pm\sqrt{6k})$
so we wind up with z=0.0896008408659, 0.749451269718 as the two fixed points. For this case, the 16 term series is accurate to about 10-11 decimal digits, and a 48 term series is accurate to 26 decimal digits. Surprisingly with a 16 term series for a=0.01, we still get ~5 decimal digits of precision.
One can also get the complex conjugate pair of fixed points for $a>-\exp(-e)$, for example for $a=0.1$ we get the following pair of complex conjugate fixed points, also accurate to about 10 decimal digits using the 16 term series above.
z=0.294596558514 +/- 0.413195411460*I
Here is a graph of the fixed points from 0.001 to 0.3. You can compare this graph to Gottfried's graph, where I have added the complex conjugate fixed points for x>exp(-e). The fixed points meet at exp(-e). We see the complex conjugate pair of fixed points to the right, and the real pair of fixed points to the left, as well as the primary real fixed point which is calculated by another method.
