$n$th root of $(a,b) \mapsto (\operatorname{gm}, \operatorname{am})$ Suppose $0 < a < b$, and let GM and AM be respectively the geometric and arithmetic means of $a$ and $b$.  Does the mapping $(a,b) \mapsto (\operatorname{GM}, \operatorname{AM})$ have a well-behaved compositional $n$th root?  In what senses might such an $n$th root be unique?
 A: Your mapping fixes the diagonal: $F(a,a)\equiv(a,a)$, so I presume that you are interested in a root $G$ that also fixes the diagonal. 
Under this natural restriction, there does not exist a twice differentiable square root $G$. 
Proof. Let $m$ be a fixed point of $F$, hence of $G$. Then $DF(m)=(DG(m))^2$ is the matrix $A$ equal to
$$\begin{pmatrix} \frac12 & \frac12 \\\\ \frac12 & \frac12 \end{pmatrix}.$$
Its only square roots are $\pm A$. But since $G$ fixes the diagonal, $DG(m)$ has an eigenvalue $1$, and therefore $DG(m)=A$.
Now, expanding at second order the identity $F(m+h)=G(G(m+h))$, and using the previous result, one obtains $AD^2G(m)=D^2F(m)$ at fixed points. But since $A$ is rank one, this implies that $h\mapsto D^2F(m)h\otimes h$ is not onto. Specifically, since $(1,-1)A=0$, one should have $D^2f(m)h\otimes h=0$ for every fixed point $m=(a,a)$ and increment $h$, where $f:=GM-AM$. But this is false, instead we have
$$D^2f(m)h\otimes h=-\frac{1}{8a}(h_2-h_1)^2.$$
QED
The proof applies to $n$th roots instead of square roots. If $G$ is an $n$th root, one still has $DG(m)=A$ at fixed point. From this, it follows that 
$$G^{(n)}(m+h)=m+Ah+\frac{1}{2}AD^2G(m)h\otimes h,$$
independently of $n$. For this calculation, it is important to notice that $m+Ah$ is a fixed point too, and that $A^2=A$. Hence the same conclusion: there does not exist a smooth $n$th root of $F$ fixing the diagonal points.
A: Here is a different (partial) answer, which does not assume any regularity of the square root $G$ of $F$ along the diagonal $a=b$. Let $M(a,b)$ denote the arithmetico-geometric mean. We have $F^{(n)}(a,b)\rightarrow(M(a,b),M(a,b))$ as $n\rightarrow+\infty$, and $M=M\circ F$. The sector $0< a< b$ is foliated by the level curves $\gamma_t=M^{-1}(t)$. By homogeneity, $\gamma_t=t\gamma_1$. Each of these curves is invariant under $F$. It is therefore natural to look for a square root $G$ that preserves every $\gamma_t$. It is enough to construct $G$ over $\gamma_1$, and then to extend it to $\gamma_t$ by homogenity: $G(ta,tb)=tG(a,b)$.
The curve $\gamma_1$ is transversal to the rays, and therefore can be parametrized by the angle $\theta\in(0,\pi/4)$. Its end point at $\pi/4$ is $(1,1)$. The restriction of $F$ over $\gamma_1$ is thus conjugated to a map $f:(0,\pi/4)\rightarrow(0,\pi/4)$. I is not hard to see that $f(\theta)> \theta$, because 
$$\frac{a+b}{2\sqrt{ab}}< \frac{b}{a}.$$
Likewise, $f'> 0$.
There remains to find a $g:(0,\pi/4)\rightarrow(0,\pi/4)$, such that $g\circ g=f$. 
A construction of $g$ could be made as follow. First find a vector field $X$ over $(0,\pi/4)$, whose flow at time $1$ is $f$. This is the hard part, for which I shall ask MO. Then set $g$ the flow of $X$ at $t=1/2$. One obtains also the $n$-th roots by taking the flow of $X$ at time $1/n$.
A: The following is more a sidenote than an answer to the problem of the fractional iteration of this map, however it introduces a connection to the half-iterate of the exponential/logarithm-map, symmetrizes the iteration and maybe interesting for that.
Define the two functions for the half-exponential and half-logarithm
$$ h(h(x)) = \exp(x)  $$
$$ g(g(x)) = \log(x) $$ 
It des not matter which base for the exponentiation we use, say we use base $ b=\sqrt 2$ which allows a real-valued solution for the half-iterate of $ b^x $ and $ log_b(x) $ (I've checked the process using that base and applying "regular fractional iteration")
Then in the original iterated map    
$$ a_{k+1} = \frac {a_k+b_k} 2 $$
$$ b_{k+1} = \sqrt {a_k b_k}  $$
substitute the initial $a_0 $ by $ A_0 = h(a_0) $ and $b_0$ by $ B_0 = h(b_0) $
Then define the iteration      
$$ A_{k+1} = g\left(\frac {h (A_k) + h (B_k)} {2} \right)  $$
$$ B_{k+1} = h(\frac {g (A_k)+g (B_k)} {2})  $$
Then $ g(A_\infty) $ and $ g(B_\infty) $ give the $\operatorname{AGM}(a,b).$
I do not see at the moment how this could be improved to allow a fractional iteration, but perhaps the idea suggests a viable direction.
A: A simple scheme (a parallel of which can be used to define the Gamma function from the factorial) is as follows: call $T$ your function, then define the functional square root as
$$T^{\frac{1}{2}}=\lim_{n\rightarrow\infty}T^{n+1}ST^{-n}$$
where $S$ is the sensible approximation of $T^{-\frac{1}{2}}$ for $b\gg a$, given by
$$S(a,b)=\Big(\frac{a^\sqrt{2}}{2^{2\sqrt{2}-5/2}b^{\sqrt{2}-1}}, \sqrt{2}b\Big)$$
Simple numerical testing shows $T^{n+1}ST^{-n}$ becoming rapidly more accurate (as $n$ grows) than $TS$ as a functional square root of $T$, and the limit above most likely converges to the only reasonable definition of $T^{\frac{1}{2}}$.
A (tedious) generalization to all fractional functional powers should be straightforward.
