What is known about the functional square root of the Riemann Zeta function? Let us consider the Riemann zeta function $\zeta(s)$ for $Re(s) > 1$: 
$$\zeta(s) := \sum_{n=1}^{\infty} \frac{1}{n^{s}} .$$
I wonder what is known about the functional square root(s) of the Riemann zeta function defined on the aforementioned domain. In other words, I'm curious about the properties of the function(s) $f$ such that $$f(f(s)) = \zeta(s). \qquad \qquad (1)$$
Questions


*

*Has a closed-form solution been found for $f$ in equation $(1)$ ? 

*If not (which I expect), have partial results been found for such a function? Properties like existence, (non)uniqueness, continuity, or results about the functional square root of the partial sums? $$f(f(s)) = \sum_{n=1}^{k} \frac{1}{n^{s}} $$

*If so, I would be grateful if you have some links to relevant articles. 

 A: One common method is to develop a series expansion about the fixed-points, that is, around where $s_\star=\zeta(s_\star)$, which occurs at $s_\star\simeq1.8338$. Now suppose that we have $s_\star=f(s_\star)$. This then let's us derive
$$\zeta'(s_\star)=f'(f(s_\star))f'(s_\star)=[f'(s_\star)]^2\\\implies f'(s_\star)=\pm\sqrt{\zeta'(s_\star)}$$
$$\zeta''(s_\star)=f''(f(s_\star))[f'(s_\star)]^2+f'(f(s_\star))f''(s_\star)=2f'(s_\star)f''(s_\star)\\\implies f''(s_\star)=\frac{\zeta''(s_\star)}{2f'(s_\star)}=\pm\frac{\zeta''(s_\star)}{2\sqrt{\zeta'(s_\star)}}$$
and so on. Since $\zeta'(s_\star)\simeq−1.374$ is negative, this gives us a non-real functional square root. This is somewhat expectable because $\zeta(s)$ behaves similarly to $s^{-1}$, which has a simple functional square root of $s^{\pm i}$.

Another simple approach is to look at rates of convergence to fixed-points. Since $\zeta$ is invertible on $(1,\infty)$, we may consider how fast $\zeta^{-n}(s)$ converges to $s_\star$. In particular, we have
$$q=\lim_{n\to\infty}\frac{\zeta^{-(n+1)}(s)-s_\star}{\zeta^{-n}(s)-s_\star}=\frac1{\zeta'(s_\star)}$$
From this, we may attempt to have
$$q^{-1/2}=\lim_{n\to\infty}\frac{\zeta^{-(n-\frac12)}(s)-s_\star}{\zeta^{-n}(s)-s_\star}=\pm\sqrt{\zeta'(s_\star)}$$
and define
$$f(s)=\lim_{n\to\infty}\zeta^n\left(s_\star+(\zeta^{-n}(s)-s_\star)q^{-1/2}\right)$$
A: In continuing @SimplyBeautifulArt's answer, I used Carleman-matrices to find a (truncated) powerseries (centered around the fixpoint) to be used for a half-iterated zeta.       
For instance, for $s_0=12$ Pari/GP gives $s_1=\zeta^{[1]}(s_0) \approx  1.00024608655 $.    
Using the Carleman-ansatz I find for    


*

*$s_{0.5}=\zeta^{[0.5]}(s_0) \approx  0.367265586625 + 0.900493834909\, î $ and 

*$s_1=\zeta^{[0.5]}(s_{0.5}) \approx  1.00024608655 + 2.87484381322 E-33 \, î $ 
where the second result is correct on more than, say, $12$ digits, which suggests that the here chosen method for determining the half-iterate is at least numerically consistent in the sense: $\zeta^{[a]}(\zeta^{[b]}(s))=\zeta^{[a+b]}(s)$.        
The whole procedure depends moreover on the implementation for the functional inverse of the zeta, written for instance as $\zeta^{[-1]}(s)$. I've no actual function for this, and used the Pari/GP-internal procedure solve() to simulate the functional inverse.    
After that, the construction, and application, of powerseries from the Carleman-matrices is in effect the Schroeder-mechanism, which uses the powerseries for the zeta centered around the fixpoint. That procedure gives by construction complex values for fractional iterates (btw. similar to the interpolation of the Fibonacci-numbers when the Binet-formula is used).

update A picture for the fractional iterates (Schroeder-method) in steps of $\Delta h=1/60$    is here:        

 
and the detail between $\zeta^{[1/2]}(12) \cdots \zeta^{[1]}(12)$ 

Remarks:     


*

*The point at iteration height $h=1$ is not exactly $1$ but about $1.00024608655$ which is not discernible in the Excel-plot.    

*The curves in the plot are Excel's interpolation-lines which are based on cubic splines      

