Dynamical System that exhibits a fold bifurcation of Limit Tori? Fold Bifurcations of a fixed points (i.e. saddle node bifurcations) and Fold bifurcations of limit cycles (i.e. when a stable limit and unstable limit cycle annihilate) are observed in plenty of different dynamical systems.
Does anyone know of a dynamical system where the equivalent happens with limit tori? An unstable and stable limit torus exist in a system, and as a parameter is tuned, these steady states approach each other until they annihilate.
 A: Yes, I think this system below, when $\lambda \in (0, 1)$, has a pair of two dimensional invariant tori, the dynamics on each of them is 'linear'.
Moreover, one of these two invariant tori exhibits hyperbolic dynamics in a tubular neighborhood,
(and it has a codimension one stable and unstable manifolds, i.e. these are three dimensional immersed submanifolds of the four dimensional phase space),
while the other torus is an asymptotically stable invariant torus. Observe the system is not just real analytic but polynomial.
When the parameter $0 < \lambda < 1$ approaches $1$, the two tori merge into one and their dynamics looks like saddle-node in transverse 2D direction. When 1 < $\lambda$
the tori disappear. For this example, assume $\omega_1$ and $\omega_2$ are some frequencies and $a = b = 0.1$
\begin{align}
&\frac{dx}{dt}\, = \, - \, \lambda \, x \, +\,\omega_1\, y \, + \,2\,x^3 \, +\,2\, xy^2 \, - \,(x^2+y^2)^2\,x \, + \, a\,\big(\, xz^2 + xw^2 - 4x\,\big)\\
&\\
&\frac{dy}{dt} \, = \, -\,\omega_1\,x \, - \, \lambda\, y \, + \,2\,x^2y \, +\,2\, y^3 \, - \,(x^2+y^2)^2\,y \, + \, a\,\big(\, yz^2 + yw^2 - 4y\,\big)\\
&\\
&\frac{dz}{dt} \, = \,\,\,\,\, 4\, z\, + \,\omega_2 \, w \, - \, z^3\, - \, zw^2 \, + \, b \, \big(\,\lambda\,z \, - \, 2\,x^2z \, - \,2\,y^2z\, + \, (x^2 + y^2)^2 \,z \,\big)\\
&\\
&\frac{dw}{dt} \, = \, -\, \omega_2\,z \, + \, 4\,w \, - \,z^2w\, - \, w^3 \, + \, b \, \big(\,\lambda\,w \, - \, 2\,x^2w \, - \,2\,y^2w\, + \, (x^2 + y^2)^2 \,w \,\big)\\
\end{align}
How can one see this? Well, the system is created so that the dynamics simplifies after performing a pair of decoupled two dimensional polar coordinate changes, i.e.
\begin{align}
&x \, = \, r \, \cos(\phi)\\
&y \, = \, r \, \sin(\phi)
\end{align}
and
\begin{align}
&z \, = \, \rho \, \cos(\psi)\\
&w \, = \, \rho \, \sin(\psi)
\end{align}
To shorten the calculations, observe that if we write the system as a pair of vector ODEs (or with complex numbers of you prefer), it is easier to perform the change of coordinates.
Let $$u = \begin{bmatrix}x \\ y\end{bmatrix} \, ,\, \,\,v = \begin{bmatrix}z \\ w \end{bmatrix}\,\, \text{ and } \,\, J = \begin{bmatrix} 0 & 1 \\ -1 & 0  \end{bmatrix} $$
as well as $|u| = \sqrt{x^2 + y^2}$ and $|v| = \sqrt{z^2 + w^2}$. Then the system above can be abbreviated as
\begin{align}
&\frac{du}{dt} \, = \, \omega_1 J \, u \, - \, \big(\,|u|^4 - 2\,|u|^2 + \lambda\,\big)\,u \, + \, a\,\big(\,|v|^2 - 4\,\big)\,u\\
&\frac{dv}{dt} \, = \,\omega_2  J \, v \, + \, \big(\,4 - |v|^2\,\big)\,v \, + \,b\,\big(\,|u|^4 - 2\,|u|^2 + \lambda\,\big)\,v
\end{align}
A straightforward calculations show that
$$|Ju| = |u| \,,\,\,\,\,\,\, |Jv| = |v| \,,\,\,\,\,\,\, u^T \,\frac{du}{dt} = \frac{d}{dt} \left( \frac{|u|^2}{2} \right)\,,\,\,\,\,\,\, v^T \,\frac{dv}{dt} = \frac{d}{dt} \left( \frac{|v|^2}{2} \right)$$
$$(J u)^T \frac{du}{dt} = |u|^2 \, \frac{d\phi}{dt} \,\,\,\,\,\,\, (Ju)^Tu = 0$$
$$(J v)^T \frac{dv}{dt} = |v|^2 \, \frac{d\psi}{dt} \,\,\,\,\,\,\, (Jv)^Tv = 0$$
so if you matrix multiply the first vector equation with $u^T$ and $(Ju)^T$ and the second equation with $v^T$ and $(Jv)^T$ the equations simplify and then it is easy to plug in the polar coordinates in each. If you prefer, (actually it may be simpler, use complex numbers). The result is the following system
\begin{align}
&\frac{dr}{dt} \, = \, - \,\big(\,r^4 - 2\,r^2 + \lambda \,\big)\,r\, + a\,\big(\,\rho^2-4\,\big)\, r \\
&\\
&\frac{d\rho}{dt} \, = \, \big(\,\rho^2-4\,\big) \, \rho \, 
- \,b \,\big(\,r^4 - 2\,r^2 + \lambda \,\big)\, \rho \\
&\\
&\frac{d\phi}{dt} \, = \, \omega_1\\
&\\
&\frac{d\psi}{dt} \, = \, \omega_2
\end{align}
Take $\lambda_1 > 0$ and $\lambda_2 > 0$ to be the two positive solutions of the polynomial $r^4 - 2r^2 + \lambda = 0$ (there are four solutions, two positive and two symmetrically negative, when $\lambda \in (0, 1)$ )
In these coordinates, the tori in question can be described as
$$T_1 = \Big\{ \, \big(r,\, \rho,\, \phi, \, \psi\big) \, \in \, \mathbb{R}^2 \times S^1 \times S^1\, \, \Big| \,\, r = \lambda_1, \,\, \rho = 2  \, \Big\}$$
$$T_2 = \Big\{ \, \big(r,\, \rho,\, \phi, \, \psi\big) \, \in \, \mathbb{R}^2 \times S^1 \times S^1\, \, \Big| \,\, r = \lambda_2, \,\, \rho = 2  \, \Big\}$$
My intuition tells me that if the ratio of the frequencies $\omega_1$ and $\omega_2$ is a Diophantine number, a highly non-rational number so to say, then the these two tori persist under small perturbations (thanks to some KAM=like theorem) and the dynamics in a neighborhood of the tori is topologically equivalent to the dynamics of this system. In other words, I think this system is structurally stable around the tori, so this phenomenon of toric saddle-node bifurcation extends to nearby polynomial or real analytic systems, close to this one.
Here is some python code that I used to plot the 2D dynamics transverse to the tori. I basically took the first two differential equations from the system in polar coordinates and plotted the phase portrait:
import numpy as np
import matplotlib.pyplot as plt


def F(x, y):
  return - x*(x**4 - 2*x**2 + lam) + a*x*(y**2-4), y*(4 - y**2) + b*y*(x**4 - 2*x**2 + lam)

r1 = 2.5
r2 = 2.5
resolution1 = 1000
resolution2 = 1000
X, Y  = np.meshgrid(np.linspace(0, r1, resolution1), np.linspace(0, r2, resolution2))

lam = 0.5
a = 0.1
b = 0.1

dXdt, dYdt = F(X, Y)

plt.streamplot(X, Y, dXdt, dYdt)
plt.axis('square')
plt.axis([0, r1, 0, r2])
plt.show()

$\lambda = 0.5$:

$\lambda = 1$:

$\lambda = 1.5$:

