One of our long-term goal is to study the origin of these satellites.  In order to do that, a model simulating the dynamical evolution under gas-drag of satellites is required.   For this purpose, we follow the approach of Beauge' and Ferraz-Mello 1993.  This it is going to be intended as a toy model, with several simplifying assumptions.  But it may be useful to understand the qualitative behavior of test particles under gas-drag.  In case of Stokes drag (see Reynolds numbers) the external force is a linear function of the relative velocity of the particle with respect to the gas (Patterson 1987) i.e.,:

Where \alpha is the ratio of gas velocity over Keplerian velocity at the same point (we assume that there is a negative radial pressure gradient, which causes the gas to have a circular velocity a bit less than the Keplerian velocity at the same point; for our purposes \alpha is 0.995), \omega is the circular angular velocity at the point (=[0, 0, sqrt(M(1)/r^3)], where r is the modulus of the position vector of the particle), and C is a nondynamical parameter which depends only on the physical properties of the gas and the density and size of the particle:
 

Where \mu is the molecular viscosity of the gas (=(m*c_m)/(3*\sigma), where m is the molecular mass, \sigma is the collisional cross section of the molecule (=2*10^-15 cm² for molecular hydrogen), and c_m is the mean thermal velocity of the gas molecule, given by sqrt((8kT)/(\pi*m)), and r_part and m_part are the particle radius and mass, respectively.  For a nebula composed solely of hydrogen molecules at a temperature of 100 K, C is given by:

We used this simple model to perform simulations on S/2000S5 and S6 over 1 Myr.   These are the results:
 
 

Note how the orbits are circularized and how, correspondingly, the value of \Theta increase with time.  For \Theta > 0.6 the satellites should get out of the Kozai resonance.  Note also the steady drop of the semimajor axis.  Results for a longer simulation (150 Myr) are shown in the following plots.