What is the link between the rotating wave approximation and the algebraic representation of a dynamical system? In analyzing a system of two coupled oscillators, I noticed a rather interesting correspondence between the so-called "rotating wave approximation" (RWA) for solving differential equations and the structure of the algebraic representation of the matrix form of those equations.
Background
The system of equations under consideration is
$$
\frac{d}{dt}
\left( \begin{array}{c} a \\ a^* \\ b \\ b^* \end{array} \right)
=
i \left(
  \begin{array}{cccc}
    - \omega_a & 0 & -g & g \\
    0 & \omega_a & -g & g \\
    -g & g & - \omega_b & 0 \\
    -g & g & 0 & \omega_b
  \end{array}
\right)
\left( \begin{array}{c} a \\ a^* \\ b \\ b^* \end{array} \right)
\tag{$\star$} \, .
$$
This system is the Hamiltonian version of the system described in a previous question, and can be derived from the Hamiltonian
$$
H = \omega_1 a^* a + \omega_2 b^* b
- g \left( a b - a^* b - a b^* + a^* b^* \right)
$$
using the appropriate Hamiltonian equations of motion
$$
\dot{a} = -i \frac{\partial H}{\partial a^*}
\quad \text{and} \quad
\dot{b} = -i \frac{\partial H}{\partial b^*}
$$
and similarly for $a^*$ and $b^*$.
The parameter $g$ characterizes the strength of coupling between the two oscillators.
Note that if $g=0$, then we have two uncoupled equations
$$
\dot a = -i \omega_a a \quad \text{and} \quad \dot b = -i \omega_b b
$$
with solutions
$$
a(t) = a(0) e^{-i \omega_a t} \quad \text{and} \quad b(t) = b(0) e^{-i \omega_b t}
$$
representing two uncoupled oscillators.
Now note that the matrix in Equation $(\star)$ can be written algebraically as
$$
  -i \sigma_z \otimes
  \left(
    g \sigma_x + \frac{d}{2} \sigma_z + \frac{s}{2} \mathbb{I}
  \right)
  - g (\sigma_y \otimes \sigma_x)
$$
where $s \equiv \omega_1 + \omega_2$ and $d \equiv \omega_1 - \omega_2$.
Rotating wave approximation
The rotating wave approximation (RWA) comes essentially from noting that while the time dependences of e.g. $a$ and $b^*$ have opposite sign and therefore lead to slow oscillation, the time dependences of e.g. $a$ and $b$ have the same sign and lead to fast oscillation.
Therefore, the terms $ab$ and $a^* b^*$ are dropped from the Hamiltonian.
Doing this leads to an equation like $(\star)$, but with a matrix given by
$$
  -i \sigma_z \otimes
  \left(
    g \sigma_x + \frac{d}{2} \sigma_z + \frac{s}{2} \mathbb{I}
  \right) \, ,
$$
i.e. the $-g (\sigma_y \otimes \sigma_x)$ term is gone.
Question
Why is the rotating wave approximation equivalent to dropping the $(\sigma_y \otimes \sigma_x)$ term in the dynamical equation matrix?
What intuitive link is there between the argument we gave for the RWA (or any other argument) and the dropping of the algebraic term $(\sigma_y \otimes \sigma_x)$?

This question has been psuedo-cross-posted to the Math site in an attempt to get a rigorous and more mathematically oriented answer.
 A: I think it may be easier to understand the system of equations if you swap the ordering of operators in ($\star$) to $(a,b,a^*,b^*)$. In this case the equation becomes
$$
\frac{d}{dt}
\left( \begin{array}{c} a \\ b \\ a^*  \\ b^* \end{array} \right)
=
-i H_{\textrm{eff}}
\left( \begin{array}{c} a  \\ b\\ a^* \\ b^* \end{array} \right)
\tag{$\star$} \, .
$$
$$
\begin{align*}
H_{\textrm{eff}} & = \sigma_z \otimes
  \left(
    g \sigma_x + \frac{d}{2} \sigma_z + \frac{s}{2} \mathbb{I}
  \right)
  - i\,  g (\sigma_y \otimes \sigma_x) \\
&=  \left(
  \begin{array}{cccc}
    \omega_1 & g & 0 & -g \\
    g & \omega_2 & -g & 0 \\
    0 & g & -\omega_1 & -g \\
    g & 0 & -g & -\omega_2
  \end{array}
\right)
\end{align*}
$$
The $ i\, g\, \sigma_y\otimes \sigma_x $ term then corresponds to the block off-diagonal part of the matrix. The first diagonal block has eigenvalues $\frac{s}{2} \pm \frac{1}{2}\sqrt{d^2 + 4 g^2}$, where $s = \omega_1 + \omega_2, d = \omega_1-\omega_2$. To get the second block's eigenvalues we send $s\rightarrow-s$. So for $s\gg |d|, g$, these blocks are split in 'energy' by approximately $s$. 
Once the block off-diagonal terms are removed, each of the two remaining blocks contains only both positive or both negative frequencies, i.e. the clockwise or counterclockwise rotating parts of each physical mode.
Therefore, one intuitive way to understand the meaning of the rotating wave approximation is that it keeps only dynamical terms that rotate in the same direction.
The 'secular approximation' for this 'Hamiltonian' corresponds to throwing away terms that are 1) small and 2) do not 'conserve energy'. In the case $s\gg |d|, g$, the term $ i\, g\, \sigma_y\otimes \sigma_x $ is small and changes the energy of a 'state' by approximately $\pm s$. This is unlike the other term proportional to $g$, which although is small does not change the energy of a 'state' by more than $\sim |g|$.
A: The general frame or these approximations is known as "secular approximation", that can be formulated as "Fast oscillating terms do not contribute to the long-term evolution of the system". Implicitly this assumes that there are non-vanishing slow terms that alone drive the long-term evolution, where "long" means at a timescale much larger than the period of oscillations. This approximation is deeply related with adiabatic aporoximation.
