Return to Answer

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.
      Orbits http://i44.tinypic.com/24pfhxy.jpg
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.

      Rotation numbers http://i41.tinypic.com/25i9bp5.jpg
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.

      Rotation numbers zoomed http://i44.tinypic.com/2hwp3yp.jpg
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).
      Mode Locked Orbits http://i40.tinypic.com/x0xx1k.jpg
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).

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).

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 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.
      Orbits http://i44.tinypic.com/24pfhxy.jpg
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.

      Rotation numbers http://i41.tinypic.com/25i9bp5.jpg
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.

      Rotation numbers zoomed http://i44.tinypic.com/2hwp3yp.jpg
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$.

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.
      Orbits http://i44.tinypic.com/24pfhxy.jpg
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.

      Rotation numbers http://i41.tinypic.com/25i9bp5.jpg
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.

      Rotation numbers zoomed http://i44.tinypic.com/2hwp3yp.jpg
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$.

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$ 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

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