My answer applies to standard continued fractions $x_0=[a_0,a_1,\dots]$ corresponding to $\alpha=1$, as well as to the full general case. As stated in my original question, the $a_n$'s are called the digits of $x_0$ in the continued fraction system (actually, continued power if $\alpha\neq 1$). The following is true for all *normal* numbers $x_0$ in that system, with the definition of normality being a generalization of the concept as we know it in base $10$ or base $2$. All numbers $x_0\in[0,1]$ but a set of Lebesgue measure zero, are normal in the continued power numeration system. Yet the set of non-normal numbers is dense.

That said, let's now answer my question. The theoretical distribution $D$ of any digit $a_n$ of any normal number $x_0$ in the power system, satisfies the following:

 1. It is unbounded, that is, $P(D=k)>0$ for any strictly positive integer $k$
 2. $P(D=k) \sim O(k^{-(1+1/\alpha)})$ as $k\rightarrow\infty$
 3. $E(D)$ is infinite if and only if $\alpha\geq 1$
 4. $Var(D)$ is infinite if and only if  $\alpha\geq \frac{1}{2}$

The correlation between successive digits thus exists if $\alpha<\frac{1}{2}$.

**Proof of the main result** (sketch)

Let $X$ be the random variable associated with the equilibrium distribution of the sequence $(x_n)$. We have

$$P(D=k)=P(m_{k+1}<X\leq m_k)=\int_{m_{k+1}}^{m_k} f_X(x)dx \sim O(m_k-m_{k+1})$$

where $m_k=(\lambda/k)^{1/\alpha}$. Since $f_X(x)$ is bounded, smooth, and away from zero in $[0,1]$, results 1 and 2 follow easily. The asymptotic result is based on the fact that $f_X$ is quite close to a uniform density on $[0,1]$. 

We also have

$$E(D)=\int_0^1 \frac{\lambda}{x^\alpha}f_X(x)dx - \int_0^1 \Big\{\frac{\lambda}{x^\alpha}\Big\}f_X(x)dx$$
where the curly brackets represent the fractional part function. The rightmost integral is finite, and the left one is infinite only if $\alpha\geq 1$. Now we proved the third result. A similar argument can be used to prove the last result. Note that for continued fractions ($\alpha=\lambda=1$) the distribution of $D$ is known exactly, and verifies all the properties 1-4 listed above.

**Functional equation for $f_X$**

Let $x'$ be a back image of $x$, in the sense that if $x_{n+1}=x$, then $x_{n}=x'$. For any $x\in [0,1]$, there are infinitely many back images $x_k'$ indexed by $k=1,2,\dots$ and denoted as 

$$x_k'=g_k(x) = \Big(\frac{\lambda}{x+k}\Big)^{1/\alpha}.$$

Now the back image of the whole interval $[x,x+\epsilon]$ (with $\epsilon$ very close to $0$) is
$$\bigcup_{k=1}^\infty \Big[g_k(x+\epsilon),g_k(x)\Big].$$
The length of the $k$-th interval in the above formula is approximately $\epsilon g'_k(x)$, where $g'$ denotes the derivative with respect to $x$. Thus, we have the following fundamental functional equation:
$$f_X(x)=\sum_{k=1}^\infty |g'_k(x)|\cdot f_X(g_k(x)) .$$

In particular, if $\alpha=\lambda=1$ (the standard continued fraction), it is easy to verify that $f_X(x)=\frac{1+x}{\log 2}$ satisfies the above formula and is thus the density associated with the equilibrium distribution. This is a well known result, and the proof provided here is probably one of the shortest ones.