$\begingroup$

Hmm, if we look at $g(x) = f(f(x)) = x^2-2$ and find the fixpoints $t_0=-1,t_1=2$ then we might assume a solution $h(x)$ for $g(x) = h(x-2)+2 $ and then from this we search for the half-iterate of $h(x)$ defining $h(x-2)+2 = w(w(x-2))+2$ and shall have a solution of the problem by $f(x)=w(x+2)-2$, hopefully having nonzero range of convergence at least in the vicinity of the fixpoint $t_1$.

With this ansatz we find that $$h(x) = g(x+2)-2 = (x+2)^2-4 = 4x+1x^2$$ could be a replacement for $g(x)$ and from this we get the leading terms for the power series of its half-iterate $$w(x) = 2 x + 1/6 x^2 - 1/90 x^3 + 1/720 x^4 - 7/32400 x^5 \\ + 161/4276800 x^6 - 391/55598400 x^7 + O(x^8) $$ by standard algebraic manipulations. (To check this insert $w(x)$ for $x$ in that leading terms)

Thus $f(x)$ can be written as a power series around (fixpoint) $2$: $$f(x) = 2+ 2 (x-2) + 1/6 (x-2)^2 - 1/90 (x-2)^3 + 1/720 (x-2)^4 - 7/32400 (x-2)^5 \\ + 161/4276800 (x-2)^6 - 391/55598400 (x-2)^7 + O(x^8) $$ (The link gives a confirmation using WolframAlpha and a $\cosh()$ /$\text{arccosh}()$-formula, see my other answer)

By the first $64$ terms it looks as if $w(x)$ has a positive range of convergence; and computing a couple of examples gives $$ \begin{array} {c|c|c|c} x &g(x)=h(1-t_1)+t_1 & f(x) = w(x-t_1)+t_1 & f(f(x))=w(w(x-t_1))+t_1\\ \hline\\ 1 & -1 & 0.17942912356439677341 & -1.0 \\ 1.5 & 0.25 & 1.0431497617870845281 & 0.25 \\ 2.5 & 4.25 & 3.0403583699367069623& 4.25 \\ \end{array}$$ while the results for $f(x)$ are only (but well) approximated because I've only the truncation to say 32 or 64 terms so far and not yet an analytical description for the coefficients which only would allow an exact result.

For the time being here is a plot taken on base of the power series for $f(x)$ truncated to 32 terms. I plotted $y=x,y=f(x),y=f(f(x))=g(x),y=g(f(x))$ to show more pattern. Unfortunately I don't know the Pari/GP-plot facility well enough to make the picture selfexplanatory. But I think one recognizes the $y=x$ straight line and the $y=g(x)$ parabola (red). The $y=f(x)$ curve (green) is that in between and the $y=g(f(x))$ (green) is that beyond the parabola:

This is additional information according to the request in the comment of @dfeuer

More Explanations

1) Generalities (Remark: $f(x)$ is meant here in general, not yet your sought function)

There is consent (see for instance L. Comtet, "Advanced Combinatorics", pg 144-148) that for power series without constant term there is a meaningful fractional iterate using the Bell-polynomials based on the formal power series for the function. One can express this more elegantly and concise in terms of Bell/Carleman matrixes (see wikipeda but for convenience I use the Carleman matrix here in transposed form) , which have a form such that for a function $f(x)$ we have $$ V(x) \cdot B = V(f(x)) \\ V(x) \cdot B^2 = V(f(f(x))) \\ \cdots\\ V(x) \cdot B^n = V(f^{[n]}(x)) \\ $$ where the notation $V(x)=[1,x,x^2,x^3,...]$ and thus $V(f(x))=[1,f(x),f(x)^2,f(x)^3,...]$ means infinite vectors of an indeterminate argument $x$.

This notation allows to isolate the required coefficients of the power series of $f(x)$ into the (infinite, lower triangular) matrix $B$ . Then for the operation of composition and self-composition (aka iteration) via the Bell-polynomials it suffices to denote powers of the matrix $B$ - because they implicitely define that required Bell-polynomial driven operations on the formal powerseries .

If the power series for $f(x) = \sum_{k=1}^\infty b_k\cdot x^k$ has no constant term, then the associated Bell-/Carleman-matrix is lower triangular, and powers of it can be computed also on the finitely truncated versions - giving exact coefficients which are valid also for the untruncated/infinite-size case.

Moreover this is also possible for fractional powers, and in particular for the square-root, which provides then the solution that you are searching.

2) How this concerns your problem (here $f(x)$ means your function)

The procedure is now to find an equivalent power series (or polynomial) for your function (let's call it $g(x)=f(f(x))$) which has no constant term - by shifting the $x$ and $y$-coordinate, such that $g(x) = x^2-2 $ which has a constant term, is rewritten as $h(x) = 4x + x^2$ which has no constant term and then we have modeled for any number of iterations $$g(x) = h(x-2) + 2 \\ g(g(x)) = h(h(x-2))+2 \\ \cdots$$.

Then for the function $h(x)$ we build the (lower triangular) Carleman-matrix $H$ which begins as $$ H= \small \begin{bmatrix} 1 & . & . & . & . & . \\ 0 & 4 & . & . & . & . \\ 0 & 1 & 16 & . & . & . \\ 0 & 0 & 8 & 64 & . & . \\ 0 & 0 & 1 & 48 & 256 & . \\ 0 & 0 & 0 & 12 & 256 & 1024 \end{bmatrix} \qquad \text{ matrix H is of infinite size} $$ and obviously the dot product with a $V(x)$-vector evaluates then to a $V(h(x))$-vector: $$ V(x) \cdot H = V(h(x))$$

Then the squareroot (let's call it $W$) of $H$ such that $W^2 = H$ can be determined by solution of an infinite equation system or by diagonalization.

One solution is again lower triangular and instead of a polynomially function like $h(x)$ it provides the coefficients of a power series. It begins like $$ W = \small \begin{bmatrix} 1 & . & . & . & . & . \\ 0 & 2 & . & . & . & . \\ 0 & 1/6 & 4 & . & . & . \\ 0 & -1/90 & 2/3 & 8 & . & . \\ 0 & 1/720 & -1/60 & 2 & 16 & . \\ 0 & -7/32400 & 1/540 & 1/30 & 16/3 & 32 \end{bmatrix} \qquad \text{ matrix W is of infinite size}$$ (You can check it by hand when you compute the square of W)

The second column gives the coefficients for the power series of $w(x)$ with the property $w(w(x-2))+2 = g(x)$ and thus gives an acceptable solution $f(x)=w(x-2)+2$.

Additional remark: Often with fractional iterates the occuring power series have very small or zero range of convergence - then additional measures must be introduced to make the solution usable and meaningful, but the powerseries for $w(x)$ here seems to have a convergence-radius of about $4$ and can be used for a not-too-small range of $x$ -values.