Let us define the following basic conditions for an iterated exponential function:
$$\exp^1(x)=e^x\tag{$\forall x$}$$
$$\exp^{a+b}(x)=\exp^a(\exp^b(x))\tag{$\forall a,b,x$}$$
I then pondered what sort of additional conditions could be applied. Using the useful inequality $e^x-1\ge x$, I considered adding the additional constraint:
$$\exp^a(x)-a\ge\exp^b(x)-b\tag{$a\ge b$}$$
which can be seen as a reasonable result of inductively applying the inequality. From this, I noticed that:
$$0=\exp^0(0)-0\le\exp^a(0)-a\le\exp^1(0)-1=0\tag{$\forall a\in[0,1]$}$$
$$\exp^a(0)=a\tag{$\forall a\in[0,1]$}$$
From this, one can define $\exp^a(0)$ for any $a$ by repeatedly using
$$\exp^{a+1}(0)=e^{\exp^a(0)}$$
One can also easily see that this implies $\exp^a(0)$ attains every real value exactly once, meaning it has a well-defined inverse. Now define the super-logarithm:
$$x=\operatorname{slog}(\exp^x(0))=\exp^{\operatorname{slog}(x)}(0)$$
and note that we can then write:
$$\exp^a(x)=\exp^{a+\operatorname{slog}(x)}(0)$$
which uniquely defines $\exp^a(x)$. That is to say, we have:
$$\exp^a(x)=\begin{cases}a,&x=0\land a\in[0,1]\\\ln(\exp^{a+1}(0)),&x=0\land a<0\\e^{\exp^{a-1}(0)},&x=0\land a>1\\\exp^{a+\operatorname{slog}(x)}(0),&x\ne0\end{cases}$$
One can then check that this satisfies the imposed inequality restriction as well as the functional equation. For the functional equation:
$$\exp^a(\exp^b(x))=\exp^{a+\operatorname{slog}(\exp^{b+\operatorname{slog}(x)}(0))}(0)=\exp^{a+b+\operatorname{slog}(x)}(0)=\exp^{a+b}(x)\tag{$x\ne0\land\exp^b(x)\ne0$}$$
The other cases are even simpler to prove. For the inequality:
$$\exp^a(0)-a=0\ge0=\exp^b(0)-b\tag{$\forall a,b\in[0,1]$}$$
For $a,b\notin[0,1]$, the result follows inductively. We can then see that
$$\exp^a(x)-a=\exp^{a+\operatorname{slog}(x)}(0)-(a+\operatorname{slog}(x))+\operatorname{slog}(x)$$
and so it follows for all $x$.
What interests me are conditions that do not seem unreasonable or meaningless that lead to similar uniqueness. And so here are my questions:
Is there a nice way to extend this to other bases? It seems the inequality for $e$ gets kind of messy if you try to extend it to other bases. And of course I'm not looking for something as trivial as "just linearly interpolate $\exp_b^a(0)$ for $a\in[0,1]$ with $\exp_b^1(x)=b^x$."
What other conditions can be imposed to produce a uniquely defined iterated exponential function (base $e$ or otherwise)?
And hopefully I didn't make any mistakes in the above definitions and proofs. $\ddot\smile$
The following answer should not spoil the satisfaction with that such a simple nontrivial criterion can be made in a meaningful-looking interpolation ansatz. Only since I've come across that "linear"-interpolation ansatz at various times, letting me unsatisfied because of the edgy nature of the resulting curve, I've recollected my thoughts to explain (and graphically display) my ansatz towards an improvement, which in the limit seems to give a completely smooth curve.
The key technique of my ansatz is to formulate the fractional iteration-part by a fractional matrix-power, and such fractional matrix-powers can be determined by diagonalization - or in the 2x2 size and $b=e$ case by matrix-logarithm. Nicely this special case gives the linear interpolation-method of the OP's ansatz! (However, only if the base to be used is $e=\exp(1)$ - see updated remark at end of post)
Of course, for the ideal, perfect case of a matrix-multiplication, that matrix should be of infinite size, representing the evaluation of a power series. Let the infinite matrix B be the Carleman-matrix associated to the function $\exp(x)$.
See the top-left of this matrix:
(The reciprocal factorials shown at the right side must be multiplied into the complete rows)
Then by the construction of B we can write $$ [1,x,x^2,x^3,...] \cdot B = [1,\exp(x),\exp(x)^2, \exp(x)^3, ...] $$ and by iteration $$ [1,x,x^2,x^3,...] \cdot B^2 = [1,\exp^{\circ 2}(x),\exp^{\circ 2}(x)^2, \exp^{\circ 2}(x)^3, ...]$$ For simplicity (and for implementing in Pari/GP for experimenting) I introduce the notation $V(x) = [1,x,x^2,x^3,...]$ in case of a finite matrix B just to the appropriate length. We have then $$ V(x) \cdot B^h = V( \exp^{\circ h}(x)$$ for each nonnegative integer $h$ . For us relevant is only that in the second column of $B^h$ there are the coefficients of the powerseries $ \exp^{\circ h}(x) = b_{h,0} + b_{h,1} x + b_{h,2} x^2 + ... $
The ideal interpolation-ansatz for $h=0..1$ would then just be the analytic solution for the fractional powers of $B$ in terms of the iteration-parameter $h$. But the attempt to diagonalization of $B$ to define the fractional iteration powerseries leads to the complex-valued Schröder-solution (by additional introduction of recentering the powerseries towards the (complex-valued) fixpoint of $t = \exp^{\circ h}(x) = \exp(x) = x $). But this is what we do not want here.
So let us avoid this ideal of fractional power of the infinite -$B$-path and try to sneak towards approximations with finite truncations through step-by-step-enlarging the smallest nontrivial truncation of $B$, let us begin with size $3 \times 3$ denoted as $B_3$ .
We look at $$ V_3(x) \cdot B_3 = [1, f_{3,1}(x) , f_{3,2}(x)] $$ A fractional $h$'th power of $B$ can now be done using diagonalization (in Pari/GP
mateigen(B)) giving first the three matrices $M,M^{-1}$ and the diagonal $D$ such that $$ B_3 = M_3 \cdot D_3 \cdot M_3^{-1} \text{ or for convenience }=M_3 \cdot D_3 \cdot W_3 $$ and then for fractional heights $0\le h \le 1$ we can determine easily $$ B_3^h = M_3 \cdot D_3^h \cdot W_3 $$ giving the fully functional form when we only use column $[,1]$ of the matrix $W_3$ $$ V(x) \cdot M_3 \cdot D_3^h \cdot W_3 [,1] = f^{\circ h}_{3,1}(x) $$The problem to define an initial interpolation-scheme for fractional $h$ and $f_{3,1}^{\circ h}(x)$ at $x=0$ is then to just evaluate this at $x=0$ and keep the result as functional expression in $h$ $$ V(0) \cdot M_3 \cdot D_3^h \cdot W_3[,1] = f^{\circ h}_{3,1}(0) $$
For instance we can do now a plot of the function
ploth(h=0,1, f(0,h)).By the term $D^h$ the iteration-height $h$ originally goes to the exponents of some coefficients of the explicated matrix-formula, but Pari/GP can immediately provide also a taylor-series solution with the usual consecutive powers of $h$. For instance for size $3 \times 3$ we get the powerseries (writing the shorter notation $g_3(0,h)$ instead of $f^{\circ h}_{3,1}(0)$) looking as $$ g_3(0,h) = 1.07602 h - 0.231565 h^2 + 0.166113 h^3 - 0.0178741 h^4 + 0.00769317 h^5 + O(h^6) $$ We see already by the first three coefficients, that for $h$ in the near of zero this is not too far from a linear function, btw.
Let's see what happens when we increase the matrix-size. $$ g_3(0,h) =1.07602 h - 0.231565 h^2 + 0.166113 h^3 - 0.0178741 h^4 + 0.00769317 h^5 + O(h^6) \\\ g_4(0,h) = 1.08989 h - 0.297210 h^2 + 0.271855 h^3 - 0.105184 h^4 + 0.0502430 h^5 + O(h^6) \\\ g_8(0,h) =1.09254 h - 0.325650 h^2 + 0.349011 h^3 - 0.225253 h^4 + 0.187433 h^5 + O(h^6) \\\ g_{16}(0,h) = 1.09188 h - 0.324852 h^2 + 0.350404 h^3 - 0.231533 h^4 + 0.201992 h^5 + O(h^6) \\\ $$ By the values of the coefficients alone it seems obvious that a determined approximation to some fixed powerseries occurs.
For the dimension $2 \times 2$ the diagonalization cannot be done (because the truncation of $B$ has two eigenvalues of same value $1$), but using the matrix-logarithm and exponentiation gives us first $ B_2^h = \text{EXP}(h \cdot \text{LOG}(B)) $ and with the coefficients in the column $[,1]$ we get $$ V(0) \cdot B_2^h [,1]= f^{\circ h}_{2,1} (0) = h = g_2(0,h) $$ which kindly is exactly your linear-interpolation ansatz.
Thus we can now compare the multitude of ansatzes starting with the linear case, going upwards to the powerseries for larger matrix-sizes as completely natural extensions of the linear case.
Here I have some plots for the family of interpolation-ansatzes, plotting $g_{s}(0,h)$ with sizes $s=2 \dots 16$ and $h=0 \cdots 1$ and by functional equation $\exp(g_s(0,h))$ and $\exp^{\circ 2} (g_s(0,h))$ .
To make the differences better visible I use $g_s(0,h)-h$ instead. The unit-interval of the linear interpolation ($s=2$) lies then on the unit-interval of the x-axis in the coordinates-system. We see then also its edgyness and the improved smothnesses of the larger-size interpolations.
Note, that for my own documentation I also inserted the term "Kneser-like" because it seems, that process of increasing the matrix-size runs towards the Kneser's solution, as I've documented elsewhere. Also in the picture I use 'dim' instead of 'size')
(the corners at $0.02$ steps are due to the resolution of my plot in that granularity)
In general the curves in the initiating interval $h=0 \dots 1$ are sinusoidal deviations from the linear shape, see this zoom
and more on the edginess a picture showing the first derivatives of the previous curves
Looking at your 1) question:
After that pictures I should mention, that the $2 \times 2$-case leads only to the "linear interpolation" if the base of exponentiation is just $e$.
If the base is different $b \gt 1$ and $b \ne e$ then the $2 \times 2$ matrix $B_{b:2}$ can be diagonalized and the generated critical interpolation function (for instance for base $b=3$) becomes $\small g_2(0,h)= 0.953713 h + 0.0448473 h^2 + 0.00140593 h^3 + 0.000033 h^4 + 0.00000062 h^5 + O(h^6)$ and is thus no more linear!
But all other properties generalize smoothly to larger bases, so I think this general ansatz is the direction you search for in question (1) at the end of your post.
At your 2. question: I've scribbled here a sort of general interpolation-scheme, valid for all bases $b>1$ which seem to converge to some "critical interpolation function" (also having a powerseries, likely with nonzero range of convergence) when the matrix-size is increased. The most obvious feature is, that the interpolations give always real-to-real-solution for real values $x$ and $h$. I moreover conjecture, that this ansatz converges towards the Kneser's solution (which has also just been designed to give a real-to-real solution) when the matrix-size grows unboundedly.
Because its general procedure uses truncated Carleman-matrixes, which are thus rather associated to polynomial functions $f_s(x)$ (the polynomial order according to matrix-size $s-1$ ) I coin the name "polynomial interpolation" so far. Don't have an immediate idea at hand, though, how to formulate some meaningful "uniqueness" for the whole packet presented here...