I have run some simple simulations of this system in Python, and it looks nothing like a nice integrable system once you zoom in on the orbits. Certainly, the orbits do not look like closed curves about the the fixed point, as suggested by the plot in the question. I'll give some numerical evidence of the chaotic behaviour of the map here, rather than actual proofs. However, I think that much of this behaviour could be demonstrated rigorously with nothing much more complicated than the the intermediate value theorem and some (very tedious) calculations with careful error bounds. Also, I think the kind of behaviour demonstrated by this map is reasonably common for area-preserving maps which are close to being integrable, although I'm not expert on this (the question just seems very interesting). When you perturb an integrable map, while retaining the area-preserving property, the orbits with irrational rotation numbers are relatively stable whereas the orbits at rational rotation numbers break up.
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;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.
Plotting the orbits of a set of initial points chosen along the x-axis demonstrates this behaviour, and the map looks very well-behaved so far.
Next, consider the line $\lbrace p_0+(x,0)\colon x\in\mathbb{R}^+\rbrace$. This should cut through the orbits, and we can plot the rotation number of $f$ on the orbit as a function of $x$. This can be done by appying $f$ some large number $n$ times, counting how many times the iterates rotate about the fixed point, and dividing by $n$. It appears to be an increasing function of $x$ going from about 0.210 to about 0.222, in agreement with the explanation above.
The rational number in this range with the smallest denominator is $3/14\approx0.2143$, which occurs at $x\approx1.119$. We can zoom into a small range about this value and calculate the rotation numbers again. I also multiply them by 14, so that it is clear where they pass through the integer value 3.
There is clear mode locking as the rotation number passes through 3/14. If true, this is inconsistent with the orbits being confined to closed curves about the fixed point! To show this, consider a point within the mode-locked region which is not a fixed point of $f^{14}$. There must be such points otherwise, by the fact that $f$ is analytic, $f^{14}$ would be constant everywhere, which is not the case (near the fixed point, for one thing). Then, take a closed ball which does not contain fixed points of $f^{14}$. Iteratively applying $f^{14}$ will coverge to a finite set of fixed points of $f^{14}$ on each orbit, which has measure zero, contradicting the area conservation property of $f$.
So, what does happen within the mode-locked regions? I chose a number of initial points in and around the mode-locked region and plotted a graph of the iterates of $f$ (rotated and rescaled to fit the graph).
The mode-locked region corresponding to rotation number 3/14 actually consists of 14 separate small regions (homeomorphic to closed balls) joined together, each of which contains a point about which iterates of $f^{14}$ rotate. These points form an orbit of period 14 for $f$, as do the points at which the small regions are joined (which look like hyperbolic periodic points).
It's not too difficult to see why this should be so. If you plot the angle through which $f^{14}$ rotates points on an orbit with rotation number less than 3/14, it will be less than 3 everywhere. As the rotation number of the orbit increases to 3/14, then the peaks of this graph increase to 3, representing the hyperbolic fixed points. Similarly, you can approach the rotation number of 3/14 from above, with the troughs the graph decreasing to the hyperbolic fixed points. These two limits give the two extreme (inner and outer) orbits of rotation number 3/14 joined at the hyperbolic fixed points. As they must be distinct orbits, they enclose the "mode-locked" region.
The region with rotation number 3/14 was easiest to check, as it has the lowest denominator, but I would expect that there is similar behaviour at other rational numbers. Furthermore, within the small mode-locked regions, where we saw that $f^{14}$ has orbits rotating about a central point, I would hazard a guess that we have similar behaviour again, giving a fractal structure to the orbits of $f$.
Let me now look at this from a theoretical point of view. It can be proven that it is not possible for the orbits of $f$ to all lie on closed curves about the fixed point. I'll first show that there are only finitely many periodic points of a given order (which is not a multiple of 3, although this constraint is maybe not necessary). This is enough to severely restrict the possible behaviour of the orbits of $f$ (the argument below should be perfectly rigorous once you fill in the bits I glossed over quickly).
Theorem: For any positive integer $n$, not a multiple of 3, there are only finitely solutions to $f^n(x)=x$.
Proof: First, let us upgrade $f$ to a 2-valued function by taking both signs for the square root, $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 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 set $S=\{x\in\mathbb{C}^2\colon x\in f^n(x)\}$ is algebraic (i.e., the zero set of a set of polynomials in x). As an affine variety, this is either of dimension 0 (finitely many values) or of positive dimension (uncountable and, in fact, unbounded). We just have to rule 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 set $M_{\pm}=(\pm1,-1;1,0)$ n times. So, 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 ${\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 $p_0\in\mathbb{R}^2$ denotes the fixed points of $f$.
Corollary: The set of closed curves about $p_0$ preserved by $f$ cannot cover all of $\mathbb{R}^2\setminus\{p_0\}$.
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
By a modification of this argument, you can show that the small orbits of $f^{14}$ in my plot must also break up. So, it is almost definite that some of the orbits of $f$ are in fact chaotic, and don't lie on any finite union of curves (I think you can prove this using the ideas above, but would take a lot more work to make it rigorous).
