Okay, I'd like to write a proposal for the solution, which depends on a single condition. I will detail the condition as we progress. But before that, I'd like to keep it vague. As to this; I'd like to argue a tad loosely; so ignore if I haven't made absolute samurai work of the \(\epsilon - \delta\) arguments. But I believe this question is answered absolutely in the affirmative. If Gottfried's construction of the sequence obeys the condition I am to release.
Note; this will be a fairly detailed post, and I want to get everything correct. As to that, it will be very god damned lengthy. I am first posting this here; and then I hope to answer the question on MO. But I don't want to answer on MO yet, as I'm still only 90% sure. This result is very deep, as to that. Expect me to describe some complex results off hand; if you have any questions please ask.
---------------------
We start by describing the set of functions \(C^1(H,H)\). This is the set of functions \(f\) such that \(f\) is holomorphic for \(\Re(z) < 0\) and sends to \(\Re(z) < 0\) (\(C^1\) is meant to be interpreted as once complex differentiable). Whereby, we call the set \(H = \{\Re(z) < 0\}\). These functions are effectively mappable to \(F : \mathbb{D} \to \mathbb{D}\), where \(\mathbb{D}\subset\mathbb{C}\) is the unit disk. And the set of functions \(C^1(H,H)\) is biholomorphic to \(C^1(\mathbb{D},\mathbb{D})\) through a linear fractional transformation. By which \(F(z) = \mu(f(\mu^{-1}(z)))\). This is important to remember, as a general construction.
The function \(\mu\) isn't too hard to find, but I'm too lazy to do it at the moment. And we've more important things to discuss.
An additional requirement, we are going to ask of our functions, is that \(f(0) = 0\), and which in the biholomorphic case \(F(1) = 1\). Such that, these functions are continuable to the edge of their domain. And the edge of that domain is a fixed point. This allows us to restrict \(\mu\). By which we have \(\mu : H \to \mathbb{D}\) while \(\mu(0) = 1\). Since we are mapping the boundary of a simply connected domain using an LFT, we are allowed to choose a single point; and then this mapping is unique.
The last real requirement we will make is that \(f\)'s maximal domain is \(H\), and that \(F\)'s maximal domain is \(\mathbb{D}\). From here, we enter into our problem. More traditionally this would be done using \(F\) and expanding divergent series (a la Ramanujan-Hardy circle method, as Tommy pointed out). We don't need to go the full analytic number theory route though. We get something stranger. But not as crazy as the infinite pochhammer; as also, Leo pointed out.
-------------------------
We begin by making a description that Knopp describes. But we do so, based on an additional criterion on our function \(F\). This is the part, I'm having trouble with; and which this is not quite a theorem, but rather a conditional theorem. But, assuming we can expand:
\[
|F(z+1) - \sum_{k=0}^N D_k z^k| \le \vartheta z^{N+1}\\
\]
So long as \(z+1 \in \mathbb{D}\)--where \(0 \le \vartheta < 1\) is not analytic, and not a proper derivative--just a constant which satisfies the asymptotic as \(z \to 0\). This DOES NOT MEAN that we can assign \(\vartheta\) as the \(N+1\)'th derivative at \(F(1)\). It means something much deeper; which I believe Gottfried's series is following. But I cannot be sure, as I haven't run the numbers myself. But assuming that Gottfried's numerical construction obeys the existence conditions to follow. We can agree.
--------------------------
The second disclaimer is made here. This will serve as the complex iteration/complex dynamics/tetration stuff so familiar in this forum. It is that (Recalling \(H = \{\Re(z) < 0\}\)):
\[
f(z) = e^z - 1\\
\]
Satisfies:
\[
f^{\circ n}(H) \to 0\,\,\text{as}\,\,n\to\infty\\
\]
This is pretty standard. Additionally we have that:
\[
f^{\circ n}(-H) \to \infty\\
\]
And for the line \(z \in i\mathbb{R}\), we are given absolute chaos. There exists a fixed point at \(0\), which we write as:
\[
f(z) = z + O(z^2)\\
\]
By then, we have an abel function \(\alpha\) about the attracting petal. We are again, following Milnor, who follows Ecalle. By which we have \(\alpha\) holomorphic for \(\Re(z) < 0\). We can hammer this home by recalling that, \(f'(z) \neq 0\). And we can always write:
\[
\alpha(f^{\circ n}(z)) - n = \alpha(z)\\
\]
By which an inverse function:
\[
f^{\circ -n} \left(\alpha^{-1}(z)\right) = \alpha^{-1}(z-n)\\
\]
And since \(f\) is non-singular \(f'(z) \neq 0\). We have a surjective map that is locally injective:
\[
\begin{align}
\alpha(z)&: H \to \mathbb{C}_{\Re(z) > 0}\\
\alpha^{-1}(z)&:\mathbb{C}_{\Re(z) > 0,|\Im(z)| < \pi} \to H\\
\end{align}
\]
Again, we need to only follow Milnor for this result. And it is fairly standard--happy to further describe this result if asked. This is, as it's known the Fatou coordinate. This is the regular iteration of \(f(z)\), in the attracting petal. Where as, if we did this in the repelling petal, we would be doing the "cheta" iteration as we recall.
But regardless, we can now get to the function we actually care about:
---------------------------
\[
g(z) = \alpha^{-1}(1/2 + \alpha(z))
\]
And this function satisfies all of the descriptions we made above. In that \(g \in C^1(H,H)\), and \(g(0) = 0\); the maximal domain of \(g\) is \(H\). This means there is a correspondent \(G : \mathbb{D} \to \mathbb{D}\). Such that \(G(1) = 1\) and it satisfies all the nice things \(g\) satisfies.
Now we know that the maximal domain of holomorphy is \(H\), by which, we have the maximal holomorphy of \(G\) is \(\mathbb{D}\). This means many tools from analytic number theory come into play. When we make guesses of series about zero; where we have a maximal domain of \(\mathbb{D}\); we are saying that the boundary of \(\mathbb{D}\) is full of singularities. So when we find a point on the boundary, that is singularity-less; (\(g(0) = 0\), \(G(1) = 1\)). We have a lot of tools to handle these exceptions.
BUT ONLY IF IT'S THE MAXIMAL DOMAIN. So what I will add again, is that \(g : H \to H\); and trying to analytically continue it to a larger domain is impossible. So, just as similarly, \(G\) is on its maximal domain at \(\mathbb{D}\). This is extraordinarily important to us, because we are summing a divergent series. And specifically for this reason. For this we can go down the rabbit hole of Ramanujan-Hardy little circle method. But we don't need to--we just need a couple of realizations from this.
Which most predominantly being, IF WE CAN FIND A SERIES:
\[
\begin{align}
g(z) &= \sum_{k=0}^\infty d_kz^k\\
G(1+z) &= \sum_{k=0}^\infty D_k z^k\\
\end{align}
\]
That they must follow similar rules. So, this isn't analytic number theory exactly. We just have to steal a couple ideas. And it's literally like \(1/100\)'th of the real analytic number theory work. Ramanujan and Hardy went really hard describing exact asymptotics of the coefficients \(D_k\)--we're just trying to pull out a simple bound. So thank god, we don't actually have to use the little circle method.
---------------------------
And here is the assumption in its complete and final form. Which, I think is perfectly reasonable; and which follows from Gottfried's coefficients. And even after reading Gottfried's construction of the coefficients--I'm sure this is valid. The trouble is I'm not 100% on constructing these coefficients on my own; so if they can be constructed in this manner--then Gottfried's solution is Borel summable. And we've constructed a new expression of \(\eta\) iteration.
But before we get there. Let's steal from analytic number theory. Let's take:
\[
G(z) = \sum_{n=0}^\infty G_n z^n\\
\]
We get the perfect estimate of:
\[
\limsup_{n\to\infty} \sqrt[n]{G_n} = 1\\
\]
Due to Cauchy's root test. Additionally, we know that:
\[
\lim_{z \to 1^-} G(z) =1\\
\]
And so, we can abel sum \(G(z)\) at \(G(1)\). Which means that:
\[
\sum_{n=0}^\infty G_n = 1\\
\]
When we Abel sum it. This is extraordinarily important because it means that \(G_n\) is Borel summable as well.
--------------------------------------------------------
Now we do the real magic, even though Abel has essentially solved our problem. We want to prove that:
\[
\sum_{n=0}^N G_n(1+z)^n = \sum_{k=0}^ND_k z^k\\
\]
And that \(D_k\) is Borel summable. The thing we can do now, because we have Abel summation--sort of has its home in "umbral calculus". But despite the shady history of "umbral calculus"--it has a strong footing in modern mathematics. But only if you can make a strong guess of the error terms in your operations. Again, which brings us to some of the ideas of analytic number theory. And here is where we invoke our "conditional theorem"--that Gottfried's series expansion satisfies the above \(\vartheta\) existence. By which. I can write:
\[
G(1+z) = \sum_{n=0}^\infty G_n(1+z)^n = \sum_{k=0}^N D_k z^k + \vartheta z^{N+1}\\
\]
For \(z\) and \(1+z\) in \(\mathbb{D}\). By which this is just \(z \in \mathbb{D} \cap (1+\mathbb{D})\). By which, if this approximation is valid, then we have our result. (I'm getting ptsd flashbacks from struggling to work through the little-circle method

)
--------------------------------------------
So let's hit it home. The values \(D_k\) can be found by abel summing \(G(1+z)\) and its derivatives. The values must satisfy this \(\vartheta\) condition. And then:
\[
\sum_{k=0}^N D_k z^k = \sum_{k=0}^N \sum_{j=0}^k \binom{k}{j} G_{k-j}z^{k-j} + \vartheta z^{N+1}\\
\]
By which we can pull out a bound of \(D_k = O(k!)\) through the binomial coefficient. This bound then applies as we conjugate back using the LFT in the beginning of this post. And we get that \(d_k = O(c^kk!)\).
I think that's as close as I'm going to get to a proof. I'd need to find something more detailed to find an exact proof. This is exactly the answer though. It'd probably take about 10 pages of strong rigour and references I was too lazy to do. But that's the answer.