## Background

In order to *really* answer this question, it is necessary to scrutinize the 'questionable manipulation' more. In particular, note that if we make a modest generalization of the power series to
$$f(a,z)=\sum_{n=0}^\infty e^{a^n} z^n$$
Then for $|z|<1$ and $|a|<1$ the series converges. Let us fix $|z|<1$ and apply your manipulation to obtain
$$F(a,z) = \sum_{k=0}^\infty \frac{1}{k!}\frac{1}{z\left(\frac{1}{z}-a^k\right)}$$
If the manipulation is correct, then for fixed $z=re^{i \theta}$ we have a pole at each $a=\sqrt[k]{r}e^{\frac{i\theta+2\pi n}{k}}$. One can numerically verify that $f$ and $F$ agree when $|a|<1,|z|<1$, thus the manipulation appears to hold inside the boundary-- and we can likely rigorously show they are equivalent by carefully dealing with the summation switch. Thus, there is a dense set of singularities on the circle $|a|=1$, so $f$ has a natural boundary. Hence we must proceed with great care since your manipulation entails evaluating a function beyond its natural boundary.

Formally, we have that
$$ a f_a(a,z)= a z f_z (a,az)$$

This is true when $|a|<1$ and continues to hold for $F$ even when $|a|>1$. Thus, $F$ at least preserves an important functional equation, and so is connected to $f$ in some way. However, if we generalize the power series to
$$g(x,a,z)=\sum_{n=0}^\infty x^{a^n} z^n$$

and apply the manipulation we obtain
$$G(x,a,z)=\sum_{k=0}^\infty \frac{1}{k!}\frac{\ln(x)^k}{z\left(\frac{1}{z}-a^k\right)}$$

It turns out that under this manipulation we end up with $g \neq G$ for some instances where $g$ is a convergent series. In particular, consider $g(x,2,-1) = \sum_{n=0}^\infty x^{2^n} (-1)^n$. In the following graph $g$ is in black and $G$ is orange.
Gottfried provides some more details on the relationship between $g$ and $G$ here: https://mathoverflow.net/a/198871/146528. The main point is that it seems that $G$ captures most of the behavior of the series, but is off by some term. Thus, I'll turn my attention to answering the question of *how* and more importantly *why* $g$ and $G$ differ in the hope that this will shed some light on how we might properly renormalize the function.

The following section includes some unrefined thoughts and theories that my study of divergent series has led me to. Thus, they should be taken with a grain of salt, since they are highly undeveloped.

## Residue at $\infty$

The residue theorem provides a link between infinite series and integration in the complex plane. For divergent series, this link is especially helpful, because we cannot evaluate the divergent series directly. Thus, if one were to establish, in full generality, a residue theorem that held even when we cannot close the contour, it would allow them to sum an extensive range of divergent series. With sufficiently strong growth conditions, we can take a tiny step toward developing a theory. A bit of this theory can be deduced from my question here: What is the relationship between $\sum_{n=0}^\infty f(n) x^n$ and $-\sum_{n=1}^\infty f(-n) x^{-n}$?. In particular, in the case that the coefficients of an alternating series are described by a complex function of exponential type less than $\pi$, we can use analytical continuation to show that the residue theorem even when we can't close the contour.

However, my studies on series with natural boundaries have illustrated that there is a piece missing in the naive theory in which we simply leave the contour unclosed and assume it contributes nothing. Consider the complex function $f(t,z) = t^{z^2}$, which is related to the coefficients of the theta function. $f$ is an entire function, yet if we integrate along the following contour (shown in grey) we obtain a non-zero result.

In fact, we can compute the value of this integral exactly by mapping from $f(z)dz \to f(\frac{1}{z})d\frac{1}{z}$, so that we are looking at $f(\frac{1}{z}) \frac{1}{z^2} dz$. Then our function looks as follows
and now computing the value of our contour comes down to integrations relating to the residue at the infinity of our function. However, the residue at infinity is zero, but our integral is non-zero. This is because the original contour didn't cover the whole plane, only the right half. The solution is to take a contour that tunnels through the essential singularity at $0$ by traveling towards it in the direction that the function is small. For instance, integrating from $-i$ up to $i$, and then connecting the two pieces by a semicircle obtains the correct result. We can connect this result directly to the series expansion of the function using a version of Cauchy's integral formula that applies to half circles. If $f(z)=\sum_{n=1}^\infty \frac{a_n}{z^n}$ then $\int_C f(z)dz = a_1\pi i-2ir\sum_{n=1}^\infty \frac{(-1)^n a_{2n}}{(2n-1)r^{2n}}$ if $C$ is a half-circle from $\theta=-\frac{\pi}{2}$ to $\theta = \frac{\pi}{2}$ with radius $r$. Thus, we may compute the residue as $\lim_{r \to 0} a_1\pi i-2ir\sum_{n=1}^\infty \frac{(-1)^n a_{2n}}{(2n-1)r^{2n}}$, which allows us to obtain that, in general, the residue picked up by the contour is $- \frac{\sqrt{\pi}}{\sqrt{-\ln(t)}}$.

Altogether, this illustrates that when we deal with functions that are large at $\infty$, we will often need to deal with extra residues at $\infty$.

## Difference $g-G$

We can do the same manipulation to obtain $G$ in a slightly different form. Using the fact that

$$\int_{-\frac{1}{2}-i\infty}^{-\frac{1}{2}+i\infty} \csc(\pi w) (za^k)^w dw = \frac{1}{1+za^k}$$
we can rewrite
$$G(x,2,-1) = \sum_{n=0}^\infty \frac{1}{k!} \frac{\ln(x)^k}{1+2^k} = \int_{-\frac{1}{2}-i\infty}^{-\frac{1}{2}+i\infty} \csc(\pi w) \sum_{k=0}^\infty \frac{\ln(x)^k}{k!}(2^k)^w = \int_{-\frac{1}{2}-i\infty}^{-\frac{1}{2}+i\infty} \csc(\pi w) x^{2^w} dw$$
If the residue theorem held this would be exactly equal to $\sum_{n=0}^\infty (-1)^n x^{2^n}$. The graph of $\csc(\pi w) x^{2^w}$ when $x<1$ sheds some light on the difference
The red contour gives $G$, and the blue contour gives $\sum_{n=0}^\infty (-1)^n x^{2^n}$. Thus, $G$ gives an incorrect answer because it picks up extra residues at infinity caused by those horizontal strips where $x^{2^n}$ is large.

In particular, integrating along the green contour picks up the largest part of the error term between $g$ and $G$, and agrees with the asymptotics derived in this answer: https://mathoverflow.net/a/198718/146528
and the strips higher up pick up the other pieces.

It is also meaningful to note that when $a<1$ the residues at infinity move from the right side of the plane to the left side, and therefore our contour doesn't pick them up, and thus the blue and red contours should be equal in this case. This agrees with the earlier computations showing $g$ and $G$ are equal when $|a|<1$ and $|z|<1$, and suggests that $G$ is indeed the correct analytic continuation when $|a|<1, |z|>1$

## Beyond the natural boundary

The difficulty in extending this function beyond its natural boundary comes from the fact that when $x>1$ the horizontal strips shift downward, so that the residue at $\infty$ cannot be seperated out from the residues we want to compute. Thus, the issue is 'how can we compute the residue at $+\infty$ without also picking up the residues at $+1, +2,+3, \dots$'. I don't know how to do this, but I'm also convinced there is some clever approach that could achieve it. I'll update this section if I manage to find something reasonable