$f(ax)=f(x)^2-1$, what is $f$?

Suppose $f(ax)=(f(x))^2-1$ and suppose that $f$ is analytic in some neighborhood of $x=0$. Expanding in power series, we get $a=1+\sqrt{5}$ or $1-\sqrt{5}$. We take positive $a$. If $f\neq{\rm const}$ then $f'(0)\neq0$ - it can be any non-zero number. After that, we can uniquely define coefficients in power series $f^{(n)}(0)/n!$ step by step using differentiation of the functional equation $$ f(0)=f(0)^2-1\ \Rightarrow\ f(0)=\frac{1\pm\sqrt{5}}2;\ 2f'(0)f(0)=af'(0)\ \Rightarrow\ f(0)=a/2;\ \ 2f''(0)f(0)+2(f'(0))^2=a^2f''(0)\ \Rightarrow\ f''(0)=\frac{2(f'(0))^2}{a^2-a};\ \ .... $$ Can $f$ be expressed in terms of some known functions?

It seems that $f$ is entire function of order $1$. Indeed, due to the Leibnitz formula, we have $$ a^nf^{(n)}(0)=\sum_{k=0}^n\binom{n}{k}f^{(n-k)}(0)f^{(k)}(0), $$ which leads to \begin{equation}\label{l1} f^{(n)}(0)=\frac{\sum_{k=1}^{n-1}\binom{n}{k}f^{(n-k)}(0)f^{(k)}(0)}{a^n-a}. \end{equation} It is true that $|f''(0)|\leq|f'(0)|^2=:C^2$, see above. Suppose that we already proved $|f^{(k)}(0)|\leq C^k$, $1\leq k\leq n-1$. Then $$ |f^{(n)}(0)|\leq \frac{C^n2^n}{(\sqrt{5}+1)^{n-1}\sqrt{5}}\leq C^n. $$ Hence, the power series converges everywhere and $|f(x)|\leq|f(0)|+e^{C|x|}-1$, $x\in\mathbb{C}$.

(The functional equation is somewhat similar to $\cos 2x=2\cos^2x-1$.

If we denote $g(x)=f(x)/f(0)$ then $$ g'(x)=g'(0)g(a^{-1}x)g(a^{-2}x)g(a^{-3}x).... $$

Also, the polynomial $$ P_n(x)=f(a^nf^{-1}(x))=P\circ...\circ P(x),\ \ \ P(x)=x^2-1 $$ is some analog of Chebyshev polynomial $T_{2^n}(x)$. The polynomials $P_n(x)$ are related to the representations of some class of infinite Lie algebras, but I forgot the story behind that...)

Remark. It is seen that $f(\lambda x)$ satisfies also the functional equation, for any $\lambda$. So, we can choose $f'(0)$ freely. Instead of this, let us choose the closest to $0$ zero as $x_0:=1$ (if zeroes exist), i.e. $f(1)=0$. Then $f(a^{-1})=1$ (not $-1$, since $x_0$ is the first zero and $f(0)>0$). Applying again the functional equation, we get $$ f(1)=0,\ f(a^{-1})=1,\ f(a^{-2})=\sqrt{2},\ ...,\ f(a^{-n})=\sqrt{1+f(a^{-n+1})},\ .... $$ Theoretically, we can recover $f(x)$ from infinite number of values $f(a^{-n})$ covergent to $f(0)=a/2$.

The next statement is wrong because, perhaps, $f$ can have complex zeros along with real zeros. I leave it just for possible improvements. Due to $f'(ax)=2a^{-1}f(x)f'(x)=(2a^{-1})^2f(x)f(a^{-1}x)f'(a^{-1}x)=...$ (see above infinite product for $g$), we can conclude that $f'(x)\neq0$ for $x\in[0,a)$ and, hence $f$ is monotonic on this interval. $f$ is also monotonic for $x\in(a,a^2)$, since $f'(ax)=2a^{-1}f(x)f'(x)$. Due to $f(a^2x)=f(x)^2(f(x)^2-2)$, we see that $x_1=a^2$ is a double zero, $x_2=a^4$ is a zero of fourth order, etc. Following the same arguments as above, there are no other zeroes. Then $$ f(x)=e^{h(x)}\prod_{n=0}^{\infty}\left(1-\frac{x}{a^{2n}}\right)^{2^n}. $$ Is $h(x)=\ln a-\ln 2$? Perhaps, something wrong can be here... but if this is true then $$ \ln f(x)=\ln a-\ln 2-\sum_{k=1}^{\infty}\frac{a^{2k}x^k}{k(a^{2k}-2)} $$ or something like that.

While the previous section was wrong, the true Hadamard expansion perhaps exist. This is still rough. Let $\{x_r\}$ be the smallest primitive zeros of $f$, such that all are other zeroes are $a^{2n}x_r$. Dentote $H(x)=\prod_r(1-x/x_r)$ - I do not know about the convergence of the product, I am trying to explain some ideas. We can follow the arguments from the wrong section above. Then, perhaps, it is true that $$ f(x)=\frac{a e^{dx}}{2}\prod_{n=0}^{\infty}H\left(\frac x{a^{2n}}\right)^{2^n} $$ with some constant $d$ (I hope it is $0$), since the order of $f$ does not exceed $1$. Substituting this identity into $f(a^2x)=f(x)^2(f(x)^2-2)$ we obtain also $$ f(x)=\sqrt{2+e^{d(a^2-2)x}(2/a)H(a^2x)}. $$ Using $f'(ax)=f'(0)\prod_{n=1}^{\infty}(2f(x/a^n)/a)$ (see above $g$), we obtain $$ f'(x)=f'(0)e^{\frac{dx}{a-1}}\prod_{n=1}^{\infty}H\left(\frac{x}{a^{2n}}\right)^{2^n-1}\prod_{n=1}^{\infty}H\left(\frac{x}{a^{2n-1}}\right)^{2^n-1}. $$ It is possible to obtain other formulas.

There is also another motivation to study $f$: $$ f(a^nx)=P\circ..\circ P(f(x)),\ \ P(x)=x^2-1. $$ Thus, we can try to analyze the stability of values $f(z)$ under the action of the group of polynomials. This question is related to dynamical systems, fractals... Maybe the dynamical systems community already studied such analytic functions?

There is an exact relation to fractals and holomorphic dynamics. Using $f(z)=P\circ...\circ P(f(a^{-n}z))$, $P(z)=z^2-1$, we obtain that if $z_0$ (located in the strip $S_n:=\{z:\ a^n\leq|z_0|<a^{n+1}\}$) is a primitive zero of $f$, i.e. $f(a^{-k}z)\neq0$, $k\geq1$, then $f(a^{-n}z_0)$ is a non-trivial root of the polynomial $P_n=P^{\circ n}$. All such non-trivia roots have the form $$ q= s_0\sqrt{1+s_1\sqrt{1+...+s_{n-1}\sqrt{1+s_n\sqrt{2}}}}, $$ where $s_j\in\{-1,1\}$. Without loss of generality, consider the case $f'(0)=1$. Then, for sufficiently large $n_0$, we have $$ \frac{a}2+a^{-n-n_0}z_0\approx f(a^{-n-n_0}z_0)= s_0\sqrt{1+s_1\sqrt{1+...+s_{n-1}\sqrt{1+s_n\sqrt{2}}}}. $$ Hence $$ a^{-n-n_0}z_0\approx s_0\sqrt{1+s_1\sqrt{1+...+s_{n-1}\sqrt{1+s_n\sqrt{2}}}}-\frac{a}2(=q) $$ and all such $q$ lying in $a^{-n_0}\leq|q|<a^{-n_0+1}$ correspond to primitive roots of $f$ lying in $S_n$. I have computed them for $n=20$ (and $n_0=4$):

It looks like a Julia set. So, zeros of $f$ (even primitive zeros) are very complex. I hope I did not make large mistakes somewhere...

It is also possible to introduce more general functions satisfying $$ f(ax)=bf(x)^2+cf(x)+d\ \ \ {\rm or\ even}\ \ \ f(ax)=Q(f(x)) $$ with some polynomial or rational $Q$. Such class of functions contain $\sin,\cos,\exp,...$

One of the main questions is still open: are there any relations between $f$ and some known functions?

Of course, any information about asymptotics, expansions, numerical results, differential equations, etc. is very welcome.

A relation with the approximation of the golden ratio. Finally, I found one relation with more or less known functions. Consider $$ g(z)=\lim_{n\to\infty}a^n\biggl(\underbrace{\sqrt{1+\sqrt{1+...\sqrt{1+z}}}}_n-\frac{a}2\biggr). $$ The function $g$ is analytic (not entire), it is considered in Paris, R. B. "An Asymptotic Approximation Connected with the Golden Number." Amer. Math. Monthly 94, 272-278, 1987. It satisfies the functional equation $$ g(z)=ag(\sqrt{1+z}). $$ Hence, $f$ is exactly the inverse function to $g$. It is useful to note that $g'$ admits an explicit representation $$ g'(z)=\frac{a}{2\sqrt{1+z}}\cdot\frac{a}{2\sqrt{1+\sqrt{1+z}}}.... $$

Edit - 04 June 2019. There are few remarks: due to the product expansion for $g'(x)$, all primitive zeros are simple and each simple zero of $f$ is primitive. Due to the comments and answers below (many thanks to the authors), the order of $f$ is $\rho=\ln 2/\ln a<1$. This value can be obtained by substituting $e^{A|z|^{\rho}}$ into $f(az)=f(z)^2-1$, which gives $a^{\rho}=2$. Because the order $\rho<1$, we have $d=0$ in the Weierstrass-Hadamard expansion. Something like this... For me, it was a bit strange to see the more or less explicit entire function $f$ whose zeros form fractals.

Edit - 05 June 2019 As mentioned above, the function $f$ has a lot of multiple roots. There is a linear transform which leaves simple roots of only. Without loss of generality, consider the case $f'(0)=1$. Recall that, see above, if $(x_r)$ are simple roots of $f$ then $f$ can be expressed in terms of $H(x)=\prod(1-x/x_r)$ as $$ f(x)=\sqrt{2+(2/a)H(a^2x)} $$ or $$ f(x)^2-1=1+(2/a)H(a^2x), $$ which leads to $$ f(x)=1+(2/a)H(ax). $$ The function $H$ has simple zeros only forming a fractal structure.

Edit - 06 June 2019 People from IMRN said that $f(az)=P(f(z))$ is a Poincare equation. They provided also somereference:

[main] P. Fatou, "Memoire sur les equations fonctionnelles", Bull. Soc. Math. Fr., 47, 161-271; 48, 33-94, 208-314 (1919).

[2] A. Eremenko and G. Levin, "Periodic points of polynomials", Ukrain. Mat. Zh., 41 (1989), 1467--1471

[3] A. Eremenko, M. Sodin, "Iterations of rational functions and the distribution of the values of Poincare functions", Teor. Funktsii Funktsional. Anal. i Prilozhen., No. 53 (1990), 18--25; translation in J. Soviet Math. 58 (1992), no. 6, 504–509

At the moment, I did not find a detailed analysis of Hadamard expansion for Poincare functions (especially for the case $P(z)=z^2-1$), but there is something in, e.g.,

[4] G. Derfel, P. Grabner, F. Vogl, "Complex asymptotics of Poincare functions and properties of Julia sets", Math. Proc. Cambridge Philos. Soc., 145 (2008), 699-718

Edit - 10 June 2019 There is a closed-form expression for $f(z)$ based on the explicit formula for inverse $g(w)=f^{-1}(w)$, see above, $$ g(w)=(w-a/2)\frac{2a}{a+2\sqrt{1+w}}\cdot\frac{2a}{a+2\sqrt{1+\sqrt{1+w}}}\dot.... $$ We can obtain explicit Vi`ete-type formulas, involving nested radicals, for all zeros of $f$. Then, Weierstrass-Hadamard factorization gives us $$ f(z)=\frac{a}{2}\prod_{\sigma\in\Sigma}\left(1+\frac{2z}{a}\prod_{n=1}^{\infty}\frac{a+2\sigma_n\sqrt{1+\sigma_{n-1}\sqrt{1+...+\sigma_1\sqrt{1}}}}{2a}\right), $$ where $$ \Sigma=\{\sigma:\mathbb{N}\to\{\pm1\},\ \ \lim_{n\to\infty}\sigma_n=1\}. $$ This factorization is one of those I was looking for.

This is more of a comment than a solution, but it is too long for a comment.

Without loss of generality you can concentrate on $a=2$. Here is a chain of substitutions starting from $$f(ax)=f(x)^2-1\rm :$$

$y=\log x$, $x=e^y$, $g(v)=f(e^v)$

This transforms $f(ax)=f(x)^2-1$ to $$g(y+\log a)=g(y)^2-1\rm .$$

$z=\frac{\log 2}{\log a}y$, $y=\frac{\log a}{\log 2}z$, $h(v)=g(\frac{\log a}{\log 2}v)$

This transforms $g(y+\log a)=g(y)^2-1$ to $$h(z+\log 2)=h(z)^2-1\rm .$$

Finally, undo the first transform:

$w=e^z$, $z=\log w$, $i(v)=h(\log v)$

This transforms $h(z+\log 2)=h(z)^2-1$ to $$i(2w)=i(w)^2-1\rm .$$

So if you find a solution to $i(2x)=i(x)^2-1$, your solution to $f(ax)=f(x)^2-1$ will be $$f(x)=i(x)^\frac{\log a}{\log 2}\rm .$$

Some partial results

If $f$ satisfies $f(a x) = f(x)^2 - 1$, then $g = f(\lambda x^n)$ satisfies $g(a^{1/n}x) = g(x^2)-1$ for every $n$ and $\lambda$. This means any solution $f$ creates a family of possible solutions. Analysis of the power series suggests that all analytic functions that satisfy $f(a x) = f(x)^2 -1$ for some $a$ are of the form $f_*(\lambda x^n)$ with $n$ an integer and have $a = a_*^{1/n}$, where ${f_*}'(0) = 1$ and $f_*(a_* x) = f_*(x)^2 - 1$. So it is reasonable to restrict to the case $f'(0) \ne 0$ provided we keep this in mind.

Next, numeric calculations. For $a = 1+\sqrt{5}$ case, the power series for $f_*$ has infinite radius of convergence and converges quickly, meaning we can calculate it numerically easily. For $x>0$, it scales roughly as $\exp(0.941x^r)$, where $r = \ln(2)/\ln(1+\sqrt{5})$ (thanks DinosaurEgg for that exponent) . For $x < 0$, it oscillates between 0 and 1 like this:

Each zero is $a^2$ times the previous and has double the order, as you found.

For $a = 1-\sqrt{5}$, the power series converges extremely slowly. Its behavior is very erf-like:

Lastly, for your infinite product representation, $h(x)$ scales as something like $0.072x^r\ln(x)$ at large positive $x$. Overall it looks like this

so I don't think it will have a simple form.

EDIT: I calculated a power series representation to 80 terms and totally forgot to check for zeros in $\mathbb C\setminus \mathbb R $. Bad news: it has them. A lot of them. For every solution to $f(z) = \pm \sqrt{2}$ and $n \ge 1$, $a^{2n}z$ is a zero. And there are infinitely many solutions to that, which means the infinite product plan probably won't work.

EDIT2: OK, so it seems the power series for the $a = 1-\sqrt{5}$ case actually does converge, just incredibly slowly--the coefficients go something like $\exp(-0.036n^{1.6})$, which grows faster than any exponential function but takes a hell of a lot of time to get there. Not sure if that exponent is $1+r$ or $(1+\sqrt{5})/2$.

$$ f(ax) = f(x)^2-1 $$

now assuming $a x > 0$ we have

$$ f(a^{\log_a(ax)})=f(a^{\log_a(x)})^2-1\Rightarrow F(u+1)=F(u)^2-1 $$

with $u = \log_a x$. Now, regarding $F(u)$, considering $F(0,y) = y$ we have using the recursion formula

$$ F(5,y) = \left(\left(\left(\left(\left(y^2-1\right)^2-1\right)^2-1\right)^2-1\right)^2-1\right)^2-1 $$

This recursion has some interesting properties.

$$ \frac{F(2,y)}{F(1,y)} = y^2 \left(y^2-2\right)\\ \frac{F(3,y)}{F(2,y)} =\left(y^2 \left(y^2-2\right)-1\right)\left(y^2 \left(y^2-2\right)+1\right)\\ \frac{F(4,y)}{F(3,y)} =y^4 \left(y^2-2\right)^2 \left(y^4 \left(y^2-2\right)^2-2\right)\\ \frac{F(5,y)}{F(4,y)} =\left(y^2-1\right)^4 \left(y^4-2 y^2-1\right)^2 \left(y^4 \left(y^2-2\right)^2 \left(y^4 \left(y^2-2\right)^2-2\right)-1\right)\\ \vdots $$

All the remainders are $-1$ and the quotients are pretty factorizable.

Also regarding the functions $F(k,y)$ they have a behavior shown in the following plot. Here are shown $F(0,y),\cdots, F(5,y)$. All curves intersect at

$$ \left(\frac{1}{2} \left(\sqrt{5}-1\right),\frac{1}{2} \left(1-\sqrt{5}\right)\right),\ \ \left(\frac{1}{2} \left(1+\sqrt{5}\right),\frac{1}{2} \left(1+\sqrt{5}\right)\right) $$

I hope this helps.