First, the function $f\colon\mathbb{R}^2\to\mathbb{R}^2$, $f(x,y)=\left(\sqrt{1+x^2}-y,x\right)$ is easily seen to be an invertible map with a single fixed point at $p_0=\left(1/\sqrt{3},1/\sqrt{3}\right)$. It has Jacobian matrix $\left(x/\sqrt{1+x^2},-1;0,1\right)$ \left(x/\sqrt{1+x^2},-1;1,0\right)$which has unit determinant, so$f$is area-preserving. At the fixed point, the Jacobian has eigenvalues$(1\pm i\sqrt{15})/4$which is not a root of unity so, near the fixed point,$f$is approximately a rotation (after a linear change of variables) by the irrational number$\theta_0=(2\pi)^{-1}\cos^{-1}(1/4)\approx0.2098$of turns. On the other hand, far from the fixed point,$f(x,y)=(\vert x\vert-y,x)+O(1)$, so$f$is approximated to leading order by$(x,y)\mapsto(\vert x\vert-y,x)$. This is integrable, with polygonal orbits (see the file linked by Sylvain Bonnot in the comments, and also the paper linked in the question, which mentions that this map has rotation number$2/9\approx0.222$). I think integrability of$f$in the limit as you go very far or very close to the fixed point is guaranteed by the fact that it is area preserving and linear on radial lines, hence reduces to a homeomorphism of the circle. 4 added 31 characters in body Proof: We can compute the rotation number of$f\vert_C$for any closed curve$C$surrounding$p_0$. For any$x\in C$, this can be calclated by counting the number of times that$f^n$rotates$x$about$p_0$, dividing by$n$, and taking the limit$n\to\infty$. This is a continuous function of$x$, which I denote by$R(x)$. As$x\to p_0$we have$R(x)\to\cos^{-1}(1/4)/(2\pi)\approx0.2098$, and$R(x)\to5/9\approx0.222$as$\Vert x\Vert\to\infty$. Choosing any rational number$p/q$(with$q$not a multiple of 3) between these limits then, by continuity of rotation numbers,$R(x)=p/q$for some$x$. As we move out radially along a line from$p_0$to infinity, one of two things can happen. (i)$R(x)=q$on a non-trivial interval. As explained above, this contradicts the area conservation property of$f$. Or, (ii)$R(x)=q$on a nowhere dense set. But, as the rotation number passes through$q$, you always get mode-locking unless$f$is conjuate to a rotation of angle$p/q$-- i.e.,$f^q(y)=y$for all$y$on the closed curve passing through$x$and preserved by$f$. However, this implies infinitely many solutions to$f^q(y)=y$. So, we can rule out both possibilities. QED 3 added proofs; edited body; added 1 characters in body I'll finally note Let me now look at this from a theoretical point of view. It can be proven that it should be is not possible for the orbits of$f$to (rigorously) determine all lie on closed curves about the fixed point. I'll first show that there are only finitely many periodic points of$f^{14}$to within a fixed given order (which is not a multiple of the computer precision3, and whether they they are hyperbolic or although this constraint is maybe not necessary). The stable and unstable manifolds This is enough to severely restrict the possible behaviour of the hyperbolic points will form smooth curves joining these points, and bound regions each containing a non-hyperbolic fixed point orbits of$f^{14}$. This would confirm f$ (the behaviour shown argument below should be perfectly rigorous once you fill in the plots drawn abovebits I glossed over quickly).
Theorem: For each unit vector any positive integer $x\in\mathbb{R}^2\setminus\{0\}$ and n$, not a multiple of 3, there are only finitely solutions to$\lambda > 0$f^n(x)=x$.
Proof: First, consider iteratively applying let us upgrade $f$ 14 times to a 2-valued function by taking both signs for the square root, $p_0+\lambda x$ f((x,y))=(\pm\sqrt{1+x^2}-y,x)$. Then,$f^n(x)$takes up to$2^n$values for any$x\in\mathbb{R}^2$, and looking at we can rewrite$f^n(x)=x$as$x\in f^n(x)$. I'll show that this actually has finitely many solutions in$\mathbb{C}^2$. Note that the number of times it rotates about set$p_0$. As S=\{x\in\mathbb{C}^2\colon x\in f^n(x)\}$ is algebraic (i.e., the zero set of a function set of $\lambda$, polynomials in x). As an affine variety, this will be a continuous function increasing from below 3 is either of dimension 0 (finitely many values) or of positive dimension (uncountable and, in fact, unbounded). We just have to above 3rule out the possibility of $S$ being unbounded. Defining the 2-valued function $g(x,y)=(\pm x-y,x)$, then $f(x)=g(x)+O(1)$. So, if there was a sequence $x_k\in S$ with $\Vert x_k\Vert\to\infty$, then $x_k/\Vert x_k\Vert\in g^n(x_k/\Vert x_k\Vert)$ up to an $O(1/\Vert x_k\Vert)$ term. Taking the limit $k\to\infty$, we have a nonzero solution to $x\in g^n(x)$. However, $g^n$ corresponds to multiplying by a matrix from the intermediate value theoremset $M_{\pm}=(\pm1,-1;1,0)$ n times. So, there for $N$ equal to one of the n-fold products $N=M_\pm M_\pm\cdots M_\pm$, we would have ${\rm det}(N-I)=0$. As these are integer matrices, we can reduce mod 2. Note that, $M_+$ and $M_-$ both reduce to $M=(1,1;1,0)$ mod 2 and that $M^3=I$ (mod 2). Therefore, $N=M^n$ which is equal to one of $M$ or $M^2$ (mod 2), as $n$ is not a multiple of 3. This implies that $\lambda(x) > 0$ so {\rm det}(M-I)=0$or${\rm det}(M^2-I)=0$(mod 2), which you can check is not the case. QED The previous theorem is enough to show that many of the curves in your plot must break up when you zoom in. Recall that$f^{14}$rotates p_0\in\mathbb{R}^2$ denotes the fixed points of $p_0+\lambda(x)x$ exactly 3 times f$. Corollary: The set of closed curves about$p_0$. Then, p_0$ preserved by $f^{14}(p_0+\lambda(x)x)=p_0+\mu(x)x$ for some f$cannot cover all of$\mu(x) > 0$. \mathbb{R}^2\setminus\{p_0\}$.
Proof: We cannot have can compute the rotation number of $\mu(x) < \lambda(x)$ f\vert_C$for all any closed curve$x\in S^1$, otherwise C$ surrounding $f^{14}$ would map the region p_0$. For any$\{p_0+\lambda x\colon x\in S^1, 0\le\lambda\le\lambda(x)\}$into its interiorC$, contradicting this can be calclated by counting the area conservation property. Similarly, number of times that $f^n$ rotates $x$ about $p_0$, dividing by $n$, and taking the limit $n\to\infty$. This is a continuous function of $x$, which I denote by $R(x)$. As $x\to p_0$ we can't have $\mu(x) > \lambda(x)$ for all R(x)\to\cos^{-1}(1/4)/(2\pi)\approx0.2098$, and$x\in S^1$R(x)\to5/9\approx0.222$ as $\Vert x\Vert\to\infty$. ThenChoosing any rational number $p/q$ between these limits then, assuming by continuity of rotation numbers, we will have $\lambda(x)=\mu(x)$ R(x)=p/q$for some$x$, in which case x$. As we move out radially along a line from $p_0+\lambda(x)x$ is p_0$to infinity, one of two things can happen. (i)$R(x)=q$on a fixed point non-trivial interval. As explained above, this contradicts the area conservation property of$x$. We should be able f$. Or, (ii) $R(x)=q$ on a nowhere dense set. But, as the rotation number passes through $q$, you always get mode-locking unless $f$ is conjuate to bound a rotation of angle $p/q$ -- i.e., $f^q(y)=y$ for all $y$ on the fixed points of closed curve passing through $f^{14}$ in x$and preserved by$f$. However, this wayimplies infinitely many solutions to$f^q(y)=y$. So, we can rule out both possibilities. ThenQED By a modification of this argument, as you can show that the Jacobian matrix small orbits of$f^{14}$is continuousin my plot must also break up. So, we can check whether these fixed points it is almost definite that some of the orbits of$f\$ are hyperbolic in fact chaotic, and don't lie on any finite union of curves (corresponding to I think you can prove this using the ideas above, but would take a trace greater than 2)lot more work to make it rigorous).