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;0,1\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][1] 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.
<br>     
![Orbits][2]
<br>
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][3] 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.
<br><br>     
![Rotation numbers][4]
<br>
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.
<br><br>     
![Rotation numbers zoomed][5]
<br>
There is clear [mode locking][6] 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).
<br>     
![Mode Locked Orbits][7]
<br>
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][8] 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$.
I'll finally note that it should be possible to (rigorously) determine all the fixed points of $f^{14}$ to within a fixed multiple of the computer precision, and whether they they are hyperbolic or not. The stable and unstable manifolds of the hyperbolic points will form smooth curves joining these points, and bound regions each containing a non-hyperbolic fixed point of $f^{14}$. This would confirm the behaviour shown in the plots drawn above.
For each unit vector $x\in\mathbb{R}^2\setminus\{0\}$ and $\lambda > 0$, consider iteratively applying $f$ 14 times to $p_0+\lambda x$ and looking at the number of times it rotates about $p_0$. As a function of $\lambda$, this will be a continuous function increasing from below 3 to above 3. So, by the intermediate value theorem, there is a $\lambda(x) > 0$ so that $f^{14}$ rotates $p_0+\lambda(x)x$ exactly 3 times about $p_0$. Then, $f^{14}(p_0+\lambda(x)x)=p_0+\mu(x)x$ for some $\mu(x) > 0$. We cannot have $\mu(x) < \lambda(x)$ for all $x\in S^1$, otherwise $f^{14}$ would map the region $\{p_0+\lambda x\colon x\in S^1, 0\le\lambda\le\lambda(x)\}$ into its interior, contradicting the area conservation property. Similarly, we can't have $\mu(x) > \lambda(x)$ for all $x\in S^1$. Then, assuming continuity, we will have $\lambda(x)=\mu(x)$ for some $x$, in which case $p_0+\lambda(x)x$ is a fixed point of $x$. We should be able to bound all the fixed points of $f^{14}$ in this way. Then, as the Jacobian matrix of $f^{14}$ is continuous, we can check whether these fixed points are hyperbolic (corresponding to a trace greater than 2).
[1]: http://www.math.harvard.edu/archive/118r_spring_05/handouts/integrable.pdf
[2]: http://i44.tinypic.com/24pfhxy.jpg
[3]: http://en.wikipedia.org/wiki/Rotation_number
[4]: http://i41.tinypic.com/25i9bp5.jpg
[5]: http://i44.tinypic.com/2hwp3yp.jpg
[6]: http://en.wikipedia.org/w/index.php?title=Circle_map&oldid=465155574#Mode_locking
[7]: http://i40.tinypic.com/x0xx1k.jpg
[8]: http://en.wikipedia.org/wiki/Hyperbolic_equilibrium_point