A more natural solution to finding the general terms of a recurrence relation in $2$ variables
There are several possible methods. A general method that I have used for similar compound means is the following. First, we want to find a recurrence for each sequence by itself. Thus, given $\,a_0 = a, \, b_0 = b, \, b_1 = \sqrt{a_1b_0}\,$ we solve for $\,a_1\,$ getting $\,a_1 = b_1^2/b_0.\,$ Similarly $\,a_2 = b_2^2/b_1.\,$ By using $\,a_2 = 2a_1b_1/(a_1+b_1)\,$ we get a quadratic in $b_2$ given $b_0$ and $b_1$ $$ b_2^2(b_0+b_1)-2 b_1^2 = 0. \tag{1} $$ Now we need an Ansatz. A little bit of numerical work suggests that the convergence of $b_n$ is linear. Thus, assume for some $\,0<r<1,\,$ without loss of generality, that $\,b_n = f(x r^n)\,$ for $\,n\ge 0\,$ where $$\,f(x) := 1 + x + c_2 x^2 + c_3 x^3 + c_4 x^4 + O(x^5). \tag{2}$$ Substitute this in equation $(1)$. The left side is $\,(r-1)(4r-1)x + O(x^2).\,$ Solving for $\,r\,$ by the hypothesis, $\,r=1/4.\,$ Solving for the other coefficients of $\,f(x)\,$ using equation $(1)$ we get $$ f(x) = 1 + x + 3/10 x^2 + 3/70 x^3 + 1/280 x^4 + 3/15400 x^5 + O(x)^6. \tag{3}$$ Looking for the denominators leads to OEIS sequence A007019
Denominators of coefficients of the Taylor series of sinh(sqrt(2*x))/(sqrt(2*x))
This gives us the solution $$ f(x) = \sinh(\sqrt{6x})/\sqrt{6x}). \tag{4}$$ Combining our results gives us $$ b_n = c\sinh(t/2^n)/(t/2^n) \tag{5}$$
Now we do something similar for $\,a_n.\,$ Thus, given $\,a_1 = 2 a_0 b_0/(a_0+b_0)\,$ we solve for $\,b_0\,$ getting $\,b_0 = a_0 a_1/(2a_0 - a_1).\,$ Similarly, $\,b_1 = a_1 a_2/(2a_1 - a_2)\,$ but since $\,b_1^2 = a_1 b_0\,$ we finally get $$ a_0(a_2^2 -4a_1a_2-4a_1^2)-a_1a_2^2 = 0. \tag{6}$$ This is a quadratic in $a_2$ given $a_0$ and $a_1$.
A little bit of numerical work suggests that the convergence of $a_n$ is linear. Thus, assume for some $\,0<r<1,\,$ without loss of generality, that $\,a_n = f(x r^n)\,$ for $\,n\ge 0\,$ where $$\,f(x) := 1 + x + c_2 x^2 + c_3 x^3 + c_4 x^4 + O(x^5). \tag{7}$$ Substitute this in equation $(6)$. The left side is $\, -(1-r)(1-4r)x + O(x^2).\,$ Solving for $\,r\,$ by the hypothesis, $\,r=1/4.\,$ Solving for the other coefficients of $\,f(x)\,$ using equation $(1)$ we get $$ f(x) = 1 + x + 6/5 x^2 + 51/35 x^3 + 62/35 x^4 + O(x^5).\tag{8}$$ The numerators are all divisible by $3$ so this suggests looking at $$ f(x/3) = 1 + 1/3 x + 2/15 x^2 + 17/315 x^3 + O(x^4). \tag{9}$$ A lookup of the numerators leads to OEIS squence A002430
Numerators in Taylor series for tan(x). Also from Taylor series for tanh(x).
which gives $$ f(x/3) = \tanh(\sqrt{x})/\sqrt{x}. \tag{10}$$ Combining our results gives us $$ a_n = c\tanh(t/2^n)/(t/2^n). \tag{11}$$
To summarize, both $\,c := \lim a_n = \lim b_n\,$ as $\,n\to\infty\,$ and $\,t\,$ depend on the initial values $\,a,b.\,$ Let $\,d := \sqrt{b^2-a^2},\,$ and $\, t := \log((b+d)/a). \,$ Then $\, c = t\ a\ b/d.\,$
Note that if $\,0<a<b\,$ then $\,t>0\,$ is real while if $\,0<b<a\,$ then $\,t\,$ is pure imaginary. Here are numerical examples for both cases:
c=1.520691, d=1.732050, t=1.316957
n a_n b_n a_n/c b_n/c (t/2^n)^2
0 1.000000 2.000000 0.657595 1.315191 1.734378
1 1.333333 1.632993 0.876793 1.073849 0.433594
2 1.468027 1.548315 0.965368 1.018165 0.108398
3 1.507103 1.527570 0.991063 1.004523 0.027099
.. ........ ........ ........ ........ ........
10 1.520691 1.520692 0.999999 1.000000 0.000001
c=1.209199, d=-1.732050, t=1.047197 i
n a_n b_n a_n/c b_n/c (t/2^n)^2
0 2.000000 1.000000 1.653987 0.826993 -1.096623
1 1.333333 1.154701 1.102658 0.954929 -0.274155
2 1.237604 1.195434 1.023459 0.988615 -0.068538
3 1.216154 1.205759 1.005751 0.997146 -0.017134
.. ........ ........ ........ ........ .........
10 1.209200 1.209199 1.000000 0.999999 -0.000001