**I'm adding a separate answer for the general question that the OP asked, which settles the question in the negative for all $n>2$ (and gives an alternate proof for $n=3$ to the one I gave above).**

Recall that the OP defined a sequence of algebraic functions $f_n$ by the rule $f_0(x) = 1$, $f_1(x) = \sqrt{x+1}$, and $f_{n+1}(x) = \sqrt{x + f_n(x)}$ for all $n\ge 1$. It was observed that $f_n$ has an elementary antiderivative for $n=0$, $1$, and $2$, and the problem was to determine whether $f_n$ has an elementary antiderivative for some $n>2$.

I am going to show that there is no elementary antiderivative of $f_n$ when $n>2$.

Assume $n>2$ (NB: This is important, because the argument below will *not* work for $n\le2$; the reader may enjoy finding where it breaks down), and let $K_n = {\mathbb C}\bigl(x,f_n(x)\bigr)$ be the elementary differential field generated by $x$ and $f_n(x)$. Then $K_n$ is the field of meromorphic functions on the normalization $\hat C_n$ of the algebraic curve $C_n$ defined by the minimal degree $y$-monic polynomial $P_n(x,y)$ that satisfies $P_n\bigl(x,f_n(x)\bigr) \equiv 0$. This minimal degree is $2^n$; for example, $P_2(x,y) = (y^2-x)^2-x-1$ and $P_3(x,y) = \bigl((y^2-x)^2-x\bigr)^2-x-1$, etc.

Since $P_{n+1}(x,y) = (P_n(x,y)+1)^2-x-1$ for $n\ge 1$ with $P_1(x,y)=y^2-x-1$, one sees, by applying the Eisenstein Criterion to $P_n(x,y)$ regarded as an element of $D[y]$ with $D$ being the integral domain ${\mathbb C}[x]$, that $P_n(x,y)$ is irreducible for all $n\ge 1$. Hence, $\hat C_n$ is connected.

It will be important in what follows to observe that $K_n$ has an involution $\iota$ that fixes $x$ and sends $f_n(x)$ to $-f_n(x)$; this is because $P_n(x,y)$ is an even polynomial in $y$. The fixed field of $\iota$ is ${\mathbb C}\bigl(x,\,f_n(x)^2\bigr)$, and the $(-1)$-eigenspace of $\iota$ is ${\mathbb C}\bigl(x,\,f_n(x)^2\bigr)f_n(x) = K_{n-1}{\cdot}f_n(x)$.

Now, the curve $C_n\subset \mathbb{CP}^2$ has only one point on the line at infinity, namely $[1,0,0]$, but the normalization $\hat C_n$ has $2^{n-1}$ points lying over this point. They can be parametrized as follows: First, establish the convention that $\sqrt{u}$ means the unique analytic function on the complex $u$-plane minus its negative axis and $0$ that satisfies $\sqrt1 = 1$ and $\bigl(\sqrt{u}\bigr)^2 = u$. Let $\epsilon = (\epsilon_1,\ldots,\epsilon_{n-1})$ be any sequence with ${\epsilon_k}^2=1$ and consider the sequence of functions $g^\epsilon_k(t)$ defined by the criteria $g^\epsilon_1(t) = \sqrt{1+t^2}$ and $g^\epsilon_{k+1}(t) = \sqrt{1+\epsilon_{n-k}t g^\epsilon_k(t)}$ for $1\le k < n$. Choose, as one may, a $\delta_n>0$ sufficiently small so that, when $t$ is complex and satisfies $|t|<\delta_n$, all of the functions $g^\epsilon_k$ are analytic when $|t|<\delta_n$. In particular, one finds an expansion
$$
g^\epsilon_n(t)
= 1+\tfrac12\epsilon_1\,t + \tfrac18(2\epsilon_1\epsilon_2-1)t^2 + O(t^3).
$$

Also, it is easy to verify that the disk in $\mathbb{CP}^2$ defined by
$$
[x,y,1] = [1,\ t g^\epsilon_n(t),\ t^2]\qquad\text{for}\quad |t|<\delta_n
$$
is a nonsingular parametrization of a branch of $C_n$ in a neighborhood of the point $[1,0,0]$. In the normalization $\hat C_n$, this is then a local parametrization of a neighborhood of a point $p_\epsilon\in \hat C_n$.
Obviously, this describes $2^{n-1}$ distinct points on $\hat C_n$.

When $x$ and $f_n$ are regarded as meromorphic functions on $\hat C_n$,
it follows that there is a unique local coordinate chart $t_\epsilon:D_\epsilon\to D(0,\delta_n)\subset \mathbb{C}$ of an open disk $D_\epsilon\subset \hat C_n$ about $p_\epsilon$ such that $t_\epsilon(p_\epsilon)=0$ and on which one
has formulae
$$
x = \frac1{{t_\epsilon}^2}
\quad\text{and}\quad
f_n(x) = \frac{g^\epsilon_n(t_\epsilon)}{t_\epsilon}
= \frac{1+\tfrac12\epsilon_1\ t_\epsilon
+\tfrac18(2\epsilon_1\epsilon_2-1)\ {t_\epsilon}^2}
{t_\epsilon}
+ O({t_\epsilon}^2).
$$
In particular, it follows that $f_n(x)$, as a meromorphic function on $\hat C_n$,
has polar divisor equal to the sum of the $p_\epsilon$ and hence has degree $2^{n-1}$. Of course, this implies that the zero divisor of $f_n(x)$ on $\hat C_n$ must be of degree $2^{n-1}$ as well.

Note that the functions $g^\epsilon_k$ satisfy $g^{-\epsilon}_k(-t) = g^{\epsilon}_k(t)$, where $-\epsilon = (-\epsilon_1,\ldots,-\epsilon_{n-1})$.
This implies that $\iota(p_\epsilon) = p_{-\epsilon}$ and that
$t_\epsilon\circ\iota = -t_{-\epsilon}$.

Now, the $2^{n-1}$ zeroes of $f_n(x)$ on $\hat C_n$ are distinct, for they are the zeros of the polynomial $q_n(x) = P_n(x,0) = (q_{n-1}+1)^2-x-1$, and the discriminant of $q_n$, being the resultant of $q_n$ and $q_n'$, is clearly an odd integer, and hence is not zero. Thus, $C_n$ is a branched double cover of $C_{n-1}$, branched exactly where $f_{n}$ has its zeros. This induces a branched cover $\pi_n:\hat C_n\to \hat C_{n-1}$ that is exactly the quotient of $\hat C_n$ by the involution $\iota$ (whose fixed points are where $f_n$ has its zeros). Since one then has the Riemann-Hurwitz formula
$$
\chi(\hat C_n) = 2\chi(\hat C_{n-1}) - B_n = 2\chi(\hat C_{n-1}) - 2^{n-1},
$$
and $\chi(\hat C_1) = \chi(\hat C_2) = 2$, induction gives $\chi(\hat C_n) =
(3{-}n)2^{n-1}$, so the genus of $\hat C_n$ is $(n{-}3) 2^{n-2} + 1$. (This won't actually be needed below, but it is interesting.)

The only poles of $x$ and $f_n(x)$ on $\hat C_n$ are the points $p_\epsilon$,
and computation using the above expansions shows that,
in a neighborhood of $p_\epsilon$, one has an expansion of the form
$$
f_n(x)\,\mathrm{d} x
- \mathrm{d}\left(f_n(x)\bigl(\tfrac12\ x + \tfrac16\ f_n(x)^2\bigr) \right)
= \left(\frac{ (1-\epsilon_1\epsilon_2) }
{4{t_\epsilon}^2}
+ O({t_\epsilon}^{-1})\right)\ \mathrm{d} t_\epsilon\ .
$$
Thus, the meromorphic differential $\eta$ on $\hat C_n$
defined by the left hand side of this equation has, at worst, double poles
at the points $p_\epsilon$ and no other poles.

Now, by Liouville's Theorem, $f_n$ has an elementary antiderivative if and only if $f_n(x)\ \mathrm{d} x$ and, hence, the form $\eta$ are expressible as finite linear combinations of exact differentials and log-exact differentials.
Thus, $f_n(x)$ has an elementary antiderivative if and only if $\eta$ is expressible in the form
$$
\eta = \mathrm{d} h + \sum_{i=1}^m c_i\,\frac{\mathrm{d} g_i}{g_i}
$$
for some $h,g_1,\cdots g_m\in K_n$ and some constants $c_1,\ldots,c_m$. Suppose that these exist. Since $\eta$ has, at worst, double poles at the $p_\epsilon$ and no other poles, it follows that $h$ must have, at worst, simple poles at the points $p_\epsilon$ and no other poles; in fact, $h$ is uniquely determined up to an additive constant because its expansion at $p_\epsilon$ in terms of $t_\epsilon$ must be of the form
$$
h = \frac{\epsilon_1\epsilon_2-1}{4t_\epsilon} + O(1).
$$

Moreover, because $\eta$ is odd with respect to $\iota$, it follows that $h$ (after adding a suitable constant if necessary) must also be odd with respect to $\iota$. This implies, in particular, that $h$ vanishes at each of the zeros of $f_n$ (which, by the argument above, are simple zeros). This implies that $h = r\,f_n$ for some $r\in K_{n-1}$ that has no poles and satisfies $r(p_\epsilon) = (\epsilon_1\epsilon_2-1)/4$ for each $\epsilon$. However, since $r$ has no poles and $\hat C_n$ is connected, it follows that $r$ is constant. Thus, it cannot take the two distinct values $0$ and $-1/2$, as the equation $r(p_\epsilon) = (\epsilon_1\epsilon_2-1)/4$ implies.

Thus, the desired $h$ does not exist, and $f_n$ cannot be integrated in elementary terms for any $n>2$.

elliptic functions, adding further evidence to the conjecture that no expression of the primitive in elementary functions can exist. $\endgroup$17more comments