Kepler orbits of settling discs
Rahul Chajwa Narayanan Menon Sriram Ramaswamy



International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru 560 089


Department of Physics, University of Massachusetts, Amherst MA 01003 USA


Department of Physics, Indian Institute of Science, Bengaluru 560 012
^*TIFR Centre for Interdisciplinary Sciences, Hyderabad 500 107 





The collective dynamics of objects moving through a viscous fluid is complex and counterintuitive <cit.>. A key to understanding the role of nontrivial particle shape in this complexity is the interaction of a pair of sedimenting spheroids  <cit.>. We report experimental results on two discs settling at negligible Reynolds number (≃ 10^-4), finding two classes of bound periodic orbits, each with transitions to scattering states. We account for these dynamics, at leading far-field order, through an effective Hamiltonian in which gravitational driving endows orientation with the properties of momentum. This leads to a precise correspondence with the Kepler problem of planetary motion for a wide range of initial conditions, and also to orbits with no Keplerian analogue. This notion of internal degrees of freedom manifesting themselves as an effective inertia is potentially a more general tool in Stokesian driven systems.

Kepler orbits of settling discs
    Sriram Ramaswamy^*
    
===============================



Particles settling in a fluid carry monopoles of force density <cit.>. In the 
Stokesian limit of Reynolds number Re → 0 they therefore manifest the hydrodynamic interaction in its strongest form. Among the consequences of this strong coupling are chaos in three-particle settling <cit.>, and the resulting statistical character of many-particle sedimentation <cit.>. Interestingly, however, the collective settling of identical spheres can be built up from two-particle processes <cit.>: i) a pair falls faster than an isolated sphere, with a horizontal drift when their separation is oblique to gravity. The reversibility of Stokes flow <cit.> ensures that the separation vector stays constant. By the same token a single apolar axisymmetric particle does not rotate, and drifts horizontally as it falls. Two non-spherical particles display a far richer sedimentation dynamics, via mutual rotation
due to a coupling between orientational and translational degrees of freedom <cit.>. In this Letter, we present experiments that classify the possible dynamical
behaviours of a settling pair of discs. We show that a symmetry-based
far-field theory, without a detailed calculation of the mutual rotation
coupling, accounts for the dynamics through the emergence of an effective
Hamiltonian for this wholly dissipative system.

Our experiments are conducted on pairs of identical discs, with radius a=0.65 cm, falling in viscous fluid (Re ∼ 10^-4) (see Methods). As shown in Figure <ref>(a) the trajectory of the centres of the discs lie in a plane. Assuming translation symmetry in this plane, the six coupled degrees of freedom can be reduced to two separation and two orientation degrees of freedom. Our observations suggest two qualitatively distinct trajectory types: scattering, in which the separation increases monotonically, and bound, in which separation and orientations oscillate with a characteristic amplitude and wavelength. The oscillatory behaviour further falls into two classes, to be discussed later. 




We ask: (i) Is there a well-defined boundary in the space of initial conditions that separates periodic and scattering (i.e. infinite-wavelength) behaviour, or do our “scattering” states simply have a wavelength longer than the container height? (ii) What determines the emergent time period and wavelength of the periodic orbits?  

Within the four-dimensional space of initial separations and orientations (Fig. <ref>), we begin with the symmetric case θ_1= 0 =θ_2, released at the same height. The resulting trajectories (Fig. <ref>) are symmetric, i.e. θ^+ ≡θ_1+θ_2=0 at all times. For small initial value x_o of the horizontal separation x, the θ_i undergo full rotations and x oscillates periodically [Fig. <ref>(a)], as noted in experiments and simulations by Jung et al. <cit.>. As x_o is increased the wavelength and amplitude of the oscillations increase sharply [Fig. <ref>(a)], until the terminal motion seems to approach the linear trajectories of isolated Stokesian discs. Finite container height makes it impossible to establish experimentally the existence of a threshold value of x_o at which the wavelength and amplitude actually diverge. A similar limitation applies to the numerical evidence for scattering orbits <cit.> using an expansion in a/R and the method of reflections <cit.>. 

Working at leading order in a/R, we  construct an effective Hamiltonian approach to the disc-settling problem and map the symmetric case to the gravitational Kepler problem, thus establishing the transition between periodic and scattering orbits. We then go on to explain the behaviours seen in asymmetric settling. We begin with an isolated settling disc: the horizontal velocity of an isolated settling disc is ẋ_̇1̇ = F αsin 2θ_1, where F is its buoyant weight and the mobility α is defined below. The tilt angle θ_1 remains constant. We can thus view the trivial evolution of x_1 and θ_1 as the Hamiltonian dynamics of a free particle with momentum θ_1 and kinetic energy proportional to cos 2θ_1. This approach also applies to the two-disc case, where θ_1,  θ_2 do not remain constant.  

For symmetric settling, retaining the lowest non-vanishing contribution in an expansion (see Supplementary Text) in a/x, ẋ = 2F αsinθ^- and θ̇^̇-̇≡θ̇_̇2̇ - θ̇_̇1̇ = 2 F γ/x^2. The proportionality constants α and γ are determined by the solution for an isolated settling spheroid <cit.>. The mobility α =- (X_A^-1 - Y_A^-1 )/12πμ a and γ =1/8πμ, where the resistance functions X_A = 8/3π and Y_A = 16/9π in the limiting case of e = √(1 - b^2/a^2)→ 1 for radius a and thickness b of the disc. The above far-field equations can be recast as Hamiltonian dynamics ẋ = ∂_θ^-ℋ,  θ̇^̇-̇ = -∂_x ℋ with   

    ℋ≡ 4F α sin^2θ^-2 + 2F γ /x
 

where 4F α sin^2θ^-2 and 2F γ / x play the roles of kinetic and potential energy respectively, with the 1/x coming from the viscous hydrodynamic kernel, not gravity. This is precisely the reduced Hamiltonian for the Kepler problem <cit.> when expressed in terms of azimuthal angle θ^- and radial coordinate x. The solution  

    1/x - 1/x_o = α/γ (cosθ^- - cosθ^-_o)
  
to the equations of motion, obtained earlier by Kim <cit.> for far-field scattering trajectories, is simply conservation of ℋ, describing both bound and scattering orbits (see Fig. <ref>c), with a transition as x_o → x_c= 4a/π. A circular Kepler orbit arises only for α=0, which is the case of a pair of identical spheres. Given the very close approach of the discs in a bound state, the far-field mapping to the Kepler problem bears up surprisingly well against experimental observations, as detailed in Fig. <ref>.



 

 A simple case of asymmetric initial conditions  consists of releasing the discs at the same height with their normal vectors perpendicular to each other, θ^+ = π/2 [Figure 3 a-c]. Once again, periodic dynamics in the orientation is observed, with the added complexity of  y oscillating between positive and negative values, and an apparent transition to unbounded orbits with increasing x_o.

The effective Hamiltonian description above provides a useful framework for understanding the dynamics resulting from a more general set of initial conditions (x_o,y_o,θ^+,θ^-_o).
A reduction to an effective two-dimensional dynamics can be achieved for asymmetric initial conditions θ^+_o ≠ 0 as well, and periodic behaviour is preserved but more complex 
(Figs. <ref> and <ref>). The resulting non-Keplerian behaviour can be understood by extending equation <ref> to incorporate the dependence of the angular velocity of the discs on the angle between the separation vector R and the external force F. To leading order in a/R, the angular velocities of discs are equal and opposite, θ̇_̇1̇ = - θ̇_̇2̇ = γ F× R /R^3. With this additional ingredient, we get the general equations of motion

    ẋ = 2 F αsinθ^-cosθ^+,       ẏ = - 2 F αsinθ^-sinθ^+
    θ̇^̇-̇ = 2 F γx/R^3,       θ̇^̇+̇ = 0.

Here x ≡ x_2 - x_1, y ≡ y_2 - y_1, θ^-≡θ_2 - θ_1 and θ^+≡θ_1 + θ_2 and α and γ are defined above (<ref>). The form (<ref>) and (<ref>) also follows on general grounds of symmetry (see Supplementary text). The conservation of θ^+ in (<ref>) constrains the dynamics of x and y to a line with slope -tanθ^+, reducing the number of variables reduces to two, thus allowing phase plane analysis. The dynamics in terms of S ≡ | R -  R_o| and θ^- (see Supplementary text) is given by Ṡ = ∂_θ^-ℋ,  θ̇^̇-̇ = - ∂_Sℋ, with effective Hamiltonian

    ℋ≡  4F α sin^2θ^-2 + 2 F γ̅(S)  R(S)

 
where γ̅(S) ≡γ( y_o - S sinθ^+)/ ( y_ocosθ^+ + x_osinθ^+) and R(S) = (S^2 + R_o^2 + 2S x_o cosθ^+ + 2S y_o sinθ^+)^1/2. Note that θ^+→ 0 yields the Keplerian limit for all initial separations (x_o,y_o).


 The Hamiltonian (<ref>) for θ^+ = π/2 implies a dynamics with y oscillating between positive and negative values, constant x, and, with increasing x_o, a transition from periodic to unbound orbits at x_c = 8 a/π (see Supplementary text). These are in accord with observations (Fig. 3 and Supplementary Videos 3 & 4), though the experiments additionally show small oscillations in x possibly arising from near-field effects and small imprecision in initial release angles.  

For both symmetric and perpendicular initial conditions, the time period diverges at the boundary between bound and scattering orbits. Assuming the ratio of thickness to radius of discs is negligibly small, we have two length scales in the problem: the radius a of the disc and the separation R between the particles. One expects the period T = (a^2μ/F) f(R/a, Re , Fr , θ^+, θ^-_o), where the scaling function f depends on the initial orientations, as well as on the Reynolds number Re = ρ U a/μ≃ 10^-4  and Froude number Fr = U/√(g a)≃ 10^-3 both of which are negligibly small. f(R_o/a) can be calculated for symmetric and perpendicular cases (see Supplementary text) in the far-field limit, whence we find that the wavelength λ∼ T F / a μ diverges more strongly (∼ A^3) for the perpendicular case than for the symmetric case (∼ A^3/2 Kepler's 3rd Law), a trend consistent with our observations (see Figure 3d). 



Rocking – a qualitatively distinct periodic behaviour analogous to libration in a pendulum, in which θ^- oscillates in a limited range – emerges for π/2 < θ^+ < π. Releasing the discs with θ_1 = π/2 and decreasing -θ_2 from π/2 (symmetric case) towards zero we experimentally capture the tumbling-rocking transition  at θ^-_o=-π/2 (see Figure 4a and 4b). Unlike in tumbling, in rocking orbits the sign of x and hence, from (<ref>), of θ̇^̇-̇, alternates as the particles interchange their relative horizontal positions. Except for the special cases of parallel and perpendicular release, rocking dynamics is best viewed in x, y and θ^- space albeit with proportional x and y displacements. Figure 4(c) shows the trajectories projected on the x-θ^- plane. The tumbling-rocking transition can once again be understood in terms of the effective Hamiltonian (<ref>) (see Supplementary text). 

Our experiments have uncovered a rich dynamics in the zero-Reynolds-number settling of a pair of identical discs, with a well-defined boundary between bound and scattering orbits and two distinct classes of periodic bound-state motion. Despite limited accuracy in locating the bound-scattering boundary, and excluding extreme situations where a disc is in the hydrodynamic shadow of another, the far-field hydrodynamic interaction offers a satisfactory and detailed understanding of the dynamics, even close to particle contact. Unexpectedly, the conservative dynamics generated by an effective Hamiltonian governs this viscosity-dominated system, with the tilt of the discs playing the role of momentum. For a large family of initial conditions the problem maps precisely to that of Kepler orbits. We also find and account for a distinct family of orbits with no planetary-orbit analogue, where the angle executes oscillations over a limited range. We expect our approach to offer insight into order and chaos in the settling dynamics of a wide class of anisotropic-particle systems.
 


METHODS

The settling experiments were carried out in a quasi-two-dimensional container with dimensions of 30 cm x 50 cm x 5 cm (Width x Height x Depth). The viscous fluid was transparent polydimethylsiloxane (silicone oil) of viscosity 60000 cSt and density 0.96 g cm^-3. The particles were made by punching out discs of diameter d=1.3 cm  from a 1mm thick aluminium sheet of density 2.7 g cm^-3. They were made smooth using sandpaper and were spray painted black. Time-lapse images were captured every 5 seconds with a Nikon D7000 DSLR using the gphoto2 commandline interface. The imported images were converted to 8-bit, and thresholded after subtracting the background. The tracking was done by fitting an ellipse to the discs, with the centroid of the ellipse giving the positions (x_i, y_i) with an error of ± 0.02  a and orientation of the major axis giving θ_i of the discs with an error of ±  0.061^∘.

Acknowledgements: This work was begun when the authors were affiliated to the TIFR Centre for Interdisciplinary Sciences, Hyderabad. SR acknowledges support from the Tata Education and Development Trust and from a J C Bose Fellowship of the SERB, India. We thank R Govindarajan for discussions. The experimental work at UMass was supported through NSF DMR-1507650.

99
url<#>1urlprefixURL 


HB1
authorJ. Happel & author H. Brenner,
titleLow Reynolds number Hydrodynamics: with special applications to particulate media, Kluwer publishers.
volume1,
  (year1983).
  
kim1
authorS. Kim & authorS.J. Karrila,
titleMicrohydrodynamics: Principles and Selected Applications, Dover publications, (year1991).

 

  bradybossis 
   authorJ.F. Brady  & authorG. Bossis,
   titleStokesian Dynamics.
   journalAnnu. Rev. Fluid Mech. volume20,
   pages111–157 (year1988)
   
   Chaikin
   authorJ.-Z. Xue et al., 
   titleDiffusion, dispersion, and settling of hard spheres, 
   journalPhys. Rev. Lett. volume69,
   pages1715-1718 (year1992).
     

andersonJ.L. Anderson,
titleColloid transport by interfacial forces.
 journalAnnu. Rev. Fluid Mech. volume21,
 pages61–99 (year1989)
 
 SR1
  authorS. Ramaswamy,
  titleIssues in the statistical mechanics of steady sedimentation.
  journalAdvances in Physics volume50,
  pages297–341 (year2001)
 
lauga
authorE. Lauga & authorT.R. Powers,
titleThe hydrodynamics of swimming
microorganisms.
journalRep. Prog. Phys.

  volume72, pages096601 (year2009).


jeffery1
  authorG.B. Jeffery,
  titleThe motion of ellipsoidal particles immersed in a viscous fluid.
  journalProc. R. Soc. Lond. volume102,
  pages161–179 (year1922)
 
 
 
chwang
 authorA.T. Chwang & authorT.Y. Wu,
 titleHydromechanics of low Reynolds number flow, Part 2.
 journalJ. Fluid Mech. volume67,
 pages787–815 (year1975)

wakiya
 authorS. Wakiya,
 titleMutual interaction of two spheroids sedimenting in a viscous fluid.
 journalJ.  Phys. Soc. Jpn. volume20,
 pages1502–1514 (year1965)

kim2
 authorS. Kim,
 titleSedimentation of two arbitrarily oriented spheroids in a viscous fluid.
 journalIntl J. Multiphase flow volume11,
 pages699–712 (year1985)
 
 shelley
 authorS. Jung, authorS.E. Spagnolie, authorK. Parikh, authorM. Shelley & authorA.K. Tornberg,
 titlePeriodic sedimentation in a Stokesian fluid.
 journalPhys. Rev. E volume74,
 pages035302 (year2006)
 
 witten1
  authorT. Goldfriend et al.,
  titleHydrodynamic Interactions between Two Forced Objects of Arbitrary Shape. I. Effect on Alignment.
  journalPhysics of Fluids volume27,
  pages123303 (year2015)
  
 Stokes1851
 authorG.G. Stokes, 
 titleOn the effect of the internal friction of fluids on the motion of pendulums.
 journalTrans. Camb. Phil. Soc.  volumeIX,
 pages8 (year1851)  
  
  janosi
   authorI.M. Jánosi et al. ,
   titleChaotic particle dynamics in viscous flows: The three-particle Stokeslet problem.
   journalPhys. Rev. E volume56-3,
   pages2858–2868 (year1997)
   
     Ekiel
     authorM.L.E. Jeżewska et al.,
     titlePeriodic sedimentation of three particles in periodic boundary conditions, 
     journalPhysics of Fluids volume17,
     pages093102-1-9 (year2005).
 
  Guazzelli
    authorÉ. Guazzelli & authorJ. Hinch, 
    titleFluctuations and Instability in Sedimentation, 
    journalAnnu. Rev. Fluid Mech. volume43,
    pages96–116 (year2011).
  
   witten2
   authorT. Goldfriend et al.,
   titleScreening, Hyperuniformity, and Instability in the Sedimentation of Irregular Objects, 
   journalPhys. Rev. Lett. volume118,
   pages158005-1-6 (year2017).
  
  
  
  
  
  
  
  Vijay
  authorK.V. Kumar, 
  titlePhD thesis, Indian Instt. of Science
  
  (year2010).
 
 jeffery2
  authorM. Stimson & authorG.B. Jeffery,
  titleThe Motion of two spheres in a viscous fluid.
  journalProc. R. Soc. Lond. volumeA111,
  pages110 (year1926)
 
crowley
 authorJ.M. Crowley,
 titleViscosity-induced instability of a one-dimensional lattice of falling spheres.
 journalJ. Fluid Mech. volume45,
 pages151–159 (year1971)
 
  LR
  authorR. Lahiri  & authorS. Ramaswamy,
  titleAre Steadily Moving Crystals Unstable?, 
  journalPhys. Rev. Lett.  volume79,
  pages1150–1153 (year1997).
 

LLmech
authorL.D. Landau & authorE.M. Lifshitz,
titleMechanics: Course of Theoretical Physics.
volume1.
 





Supplementary text




§ ERROR IN ELLIPSE FITTING

To find the error associated with fitting an ellipse in ImageJ, we track a single settling disc in the scattered state, where it follows a straight line path as shown in Figure 5 (a) and (b) . In Figure 5(c) residual with respect to linear fitting gives a measure of error. The measured root mean squared error is 0.012cm or 0.019 a in terms of the radius a of disc. Similarly Figure 5(d) shows the linear fit of angle of the major axis of the fitted ellipse with time giving a slope of -0.0006 and its residual in Figure 5(e) gives the root mean squared error of  0.061 degrees.  










§ LOG-LOG PLOT FOR SYMMETRIC DIVERGENCE

To determine the critical value of initial separation d for which the wavelength λ, amplitude A and time period T diverges we plot log of λ,A and T as a function of log(x_c - d/0.65) with x_c being a fitting parameter which simultaneously minimizes the difference between slope of individual curves. This minimisation gives x_c = 1 which is used for this plot. 






   


§ USING SYMMETRIES


The system as shown is figure <ref> has three dynamical variables:

	
  * K⃗^1 , K⃗^2 , axisymmetric axis of the particles 1 and 2. 
	
  * R⃗^12 separation vector pointing form particle 1 to 2. 

The equation of motion for the pair of particles is:

    V_i = M_ijF^j

where the left side is the time rate of change of the dynamical variables, M_ij is the mobility tensor and F^j is the external driving force directed vertically downwards. We construct a simple form of the second order tensor M_ij allowed by the following symmetry of this system:

	
  * Time reversal symmetry, i.e. reversal of both V and F.
	
  * Apolarity of particles, K⃗⟶ - K⃗, transform any one or both the particles.
	
  * Particle exchange symmetry
	
  * Rotational symmetry in the perpendicular subspace, x ⟶ -x , K_x⟶ -K_x and R_x^12⟶ -R_x^12.

Given the above symmetries the equation of motion for R_i looks like

    d R_i^12/ dt = M^T_ijF^j

The simplest form of M^T_ij which satisfies the above symmetries is

    M_ij = α K_i^1K_j^1  +  β K_i^2K_j^2

where the coefficients α and β can be arbitrary even scalar functions of K_lR_l. From the condition R⃗^12 = - R⃗^21, we get the relation β = -α. The spatial part of the dynamics is

    d R_x^12/dt = |F| α_x(K_x^1K_y^1  -   K_x^2K_y^2).

 Substituting K_x^i = sinθ_i and K_y^i = cosθ_i, where i = 1,2 is the particle index, gives

    d R_x^12/dt    = |F| α_x/2(sin 2θ_1  -  sin 2θ_2) = |F|α_x sin(θ_1 - θ_2)cos(θ_1 + θ_2).

and similarly for the y component,

    d R_y^12/dt    = |F| α_y(cos^2θ_1  -  cos^2θ_2) = |F|α_y sin(θ_2 - θ_1)sin(θ_1 + θ_2).

 Since the angular velocity is directed perpendicular to the plane of the dynamics, the simplest equation for the orientation angle allowed by the symmetries is

    (d θ_1/ dt )_i = γ ϵ_ijk R_j^12 F^k.

Since R^21 = - R^12 , we get 
    d θ_1/ dt = - d θ_2/ dt

This condition is not true in general and breaks down for non-symmetric configurations if we allow rotational mobility to depend on K⃗. For example, term like ϵ_ilmK_lR_mR_j K_kR_k is allowed by symmetries but is not odd under R⃗⟶ -R⃗. Also note that, the above prescription does not give us the exact form of the coefficients α_x, α_y and γ and determining their value requires hydrodynamic calculation.



§ DYNAMICS FROM A FAR-FIELD ANALYSIS

The leading order translational response of a spheroid to an external force F is

    U_1=  [X_A^-1 K K + Y_A^-1(δ -  K K)]  F/6πη a

Here K is the orientation vector, X^A and Y^A are resistance functions for spheroids <cit.>. The second spheroid generates a velocity field which is calculated using the distribution of singularities at r_2. For prolate case, it is a line distribution between the focal points:

    v_2(r) = F_2·∫_-k_2^k_2{ 1   +   (k_2^2 - r_2^2)(1-e_2^2)^2/4e_2^2∇^2}𝒢(r -r_2)/8πη  dr_2

Here 𝒢(r)_ij/8πη is the Green's function for Stokes flow and k = a e with semi-major axis a and eccentricity e. For the case of oblate spheroid the line integral is carried along a complex focal length (k ⟶ ik), which is equivalent to a distribution of singularities on a disc of radius k in the perpendicular real subspace <cit.>. Vorticity of v_2 evaluated at the singularity distribution of the first particle, placed at r_1, gives the leading contribution to the rotation of the first particle:

    ω_1 = 1/8πη∫_-k_1^k_13/4k_1^3  dr_1 ∫_-k_2^k_2  dr_2/2k_2 (k_1^2 - r_1^2)  F_2× r_12/|r_12|^3

here, r_12 = r_1 - r_2. When particles are far apart we can take the integrand F_2× r_12 / r_12^3 outside the integral and the remaining integral is unity. If R is pointing from the centroid of second particle to the centroid of first, the far-field approximation gives

    ω_1 =  F_2× R/8πη R^3



Equations of motion



Let the degrees of freedom of i^th particle be (x_i,y_i,θ_i), where i=1,2. This gives six coupled equations of motion:

    ẋ_̇1̇ = Fαsin 2θ_1      ẋ_̇2̇ = Fαsin 2θ_2
    ẏ_̇1̇ = Fβcos^2θ_1      ẏ_̇2̇ = Fβcos^2θ_2
    θ̇_̇1̇ = Fγ(x_1 - x_2)/R^3      θ̇_̇2̇ = Fγ(x_2 - x_1)/R^3


Here, α = - 1/12π a [1/X_A - 1/Y_A] , β = - 1/6π a [1/X_A - 1/Y_A] = 2α and γ = 1/8π



 Note that we have written these equations in the settling frame of the particles. In the relative coordinates: x ≡ x_2-x_1 , y ≡ y_2 - y_1, θ^-≡θ_2 - θ_1 and θ^+≡θ_2 + θ_1: 

    ẋ = 2 F αsinθ^-cosθ^+      ẏ = - 2 F αsinθ^-sinθ^+      θ̇^̇-̇ = 2 F γx/R^3      θ̇^̇+̇ = 0



§ EFFECTIVE HAMILTONIAN

From equation <ref> we get the relation between x and y

    d x/d y = - 2α/βθ^+ = -θ^+;       
    x - x_o = -θ ^+ (y - y_o).

Here, x_o and y_o are the initial x and y respectively. The conservation of θ^+ can be used to rewrite equations <ref> in terms of the arc length S along the line in the x-y plane (see figure <ref>) and θ^-.


    S = |R⃗- R⃗_⃗o⃗| = √((x - x_o)^2 + (y - y_o)^2) = (x - x_o) θ^+


    d S/d t = θ^+d x/d t = 2Fαsinθ^- =  ∂/∂θ^-(-  2F αcosθ^-)

Using equation <ref> and the constraint imposed by equation <ref> we get  

    d θ^-/d t    =  2 F γ x/(x^2 + {y_o - tanθ^+(x - x_o) }^2)^3/2  = -2 F γd /d xy_o - tanθ^+(x-x_o)/(y_o + x_otanθ^+) R

Writing <ref> in terms of S and θ^- gives

    d θ^-/d t = -2 F γθ^+d /d Sy_o - tanθ^+(x-x_o)/(y_o + x_otanθ^+)|S⃗ + R⃗_⃗o⃗| = - d /d S2Fγ(y_o - S sinθ^+)/(y_ocosθ^+ + x_osinθ^+)|S⃗ + R⃗_⃗o⃗|




§ PERPENDICULAR SETTLING


Substituting the initial horizontal separation y_o = 0 and θ^+ = π/2 in equation <ref>, <ref>, <ref> and solving for x and y in terms of θ^- gives

    x = x_o.

and

    γ x_o  dy/(x_o^2 + y^2)^3/2 = -α sinθ^-  dθ^-

which upon integration and substituting the value α/γ = π/8a for disc, gives

    y = ±cosθ^-  x_o/√((8a/π x_o)^2 - cos^2θ^-).

This simplification when θ^+ = π/2 can be made use of for a wide range of initial conditions.



§ PERIOD OF THE ORBITS

(i) θ^+→ 0



For θ^+→ 0 with arbitrary initial separation R_o = √(x_o^2 + y_o^2) and initial angle θ^-_o, the resulting trajectory becomes ellipses if separation R and angle difference θ^- is mapped to the radial coordinate and azimuthal angle respectively. 

    l/R =  1 + ecosθ^-


where the latus rectum 2l, eccentricity of the orbit e, semi-major axis p and semi-minor axis q are

    l = 8a/π/8a/ π R_o - cosθ^-_o ;        e = 1/8a/π R_o -cosθ^-_o ;        p = R_max R_min( 1/R_o - πcosθ^-_o/8a);        q = (R_max R_min)^1/2 .

here R_max = 8a/π/8a/π R_o -  cosθ^-_o -1 is the maximum separation and R_max = 8a/π/8a/π R_o -  cosθ^-_o + 1 is the minimum separation for a given initial separation R_o and angle θ^-_o. 


  §.§.§ Symmetric

The trajectories becomes mirror symmetric when y_o = 0 making R = x at all times. From equation <ref> we get θ̇^̇-̇ = 2Fγ/x^2 which is equivalent to the angular momentum conservation R^2θ̇^̇-̇ = 2Fγ. Using this conservation we can find the time period T of the orbit by an integration which amounts to calculating area of the ellipse π p q.

    T = 8 π^2μ/FR_max^3/2R_min^3/2(π/8a + 1/x_o - 1/x_c)

here x_c = 8a/π(cosθ^-_o + 1) is the critical value of initial x_o at which the amplitude R_max diverges. Time period can be written in the following form

    T = a^2μ/F f(x_o/a, θ^-_o)

where the scaling function

    f(x_o/a, θ^-_o) = 8^3( 8a/π x_o - cosθ^-_o) /( 8a/π x_o - cosθ^-_o -1 )^3/2( 8a/π x_o - cosθ^-_o + 1 )^3/2

When x_o approaches x_c the time period diverges with amplitude as T ∼ A^3/2. Since wavelength λ scales as λ∼ T F / a μ it diverges the same way as T.



(ii) θ^+→π/2



The solution of equation <ref> and <ref> in the limit θ^+→π/2 is

    y_o-y/R = π x_o/8a (cosθ^- - cosθ^-_o)



  §.§.§ Perpendicular

when the initial angle θ^-_o = π/2 and initial y_o =0 the trajectories takes the form of ellipse: l/ r =  1 + ecosϕ if we map the square of separation R^2 to the radial coordinate r  and 2 θ^- - π to the  azimuthal angle ϕ. Here the latus rectum 2l, eccentricity of the orbit e, semi-major axis p and semi-minor axis q are

    l = 2 a^2/2 a^2/ x_o^2 - (π/8)^2 ;        e = (π/8)^2/2 a^2/ x_o^2 - (π/8)^2 ;        p = R_max^2/2( 2(8a/π x_o)^2 - 1 );        q = 8a/πR_max .

here R_max = x_o/√((8a/π x_o)^2 -1) is the maximum separation or amplitude of the oscillation. From equation <ref> we have

    R^3dθ^-/dt = 2Fγ x_o

rewriting it in terms of r and ϕ

    1/2√(r)r^2dθ^- = 2Fγ x_odt

integrating the left hand side by parts

    1/√(r)∫r^2/2dϕ    -∫√(r)/8(π/8a)^2 sinϕ  dϕ ∫r^2/2dϕ  =   2Fγ x_o T

over one complete cycle ϕ goes from 0 to 2π . The indefinite integral over r^2dϕ/2 is the area of sector of an ellipse which does not simplify <ref> any further. To get a closed form of period as a function of x_o we approximate the area of an eliptical sector by its supremum value p   q ϕ/2 . The second integral on the left hand side can be written completely in terms of ϕ using equation <ref>

    π  p   q/x_o    - p q/16 ∫_0^2π(π/8a)  ϕ sinϕ/√((8a/π x_o)^2 - 1/2 + cosϕ/2)  dϕ   =   2Fγ x_o T

integrating the left hand side again by parts gives

    π  p   q/x_o    - π p q/32 a ∫_0^2π  dϕ  √((8a/π x_o)^2 - 1/2 + cosϕ/2)  =   2Fγ x_o T

the above integral has the form of complete elliptic integral E[ϕ/2 , (π x_o/8a)^2]. This can be solved for time period T in terms of x_o by inseting the value for p and q

    T = 32  a x_o μ/F  f(x_o/a)

here the scaling function f(x_o/a) is

    f(x_o/a) =  2 (x_c/x_o)^2 -1 /[ (x_c/x_o)^2 - 1 ]^3/2 ( π  -  (x_o/x_c)^2  E[(x_o/x_c)^2] )

The critical value of x_o at which f(x_o/a) diverges is x_c = 8a/π. As x_o→ x_c the time period T diverges as

    T = 256  a^2 μ/F[ (x_c/x_o)^2 - 1 ]^3/2

in terms of amplitude A, time period goes as T ∼ A^3 as x_o→ x_c.



§ ROCKING TRAJECTORIES

For initial θ^-_o<π/2 both the disc rotates in a range of angles which defines the amplitude of angular oscillation. figure <ref> shows the trajectories of one of the disc for θ^-_o<π/2 which exhibits oscillation in a range of angles (red) and θ^-_o>π/2 showing tumbling orbits (blue).

Using equation <ref> we construct the dynamics in x-θ^- plane:

    x dx/(x^2 + y^2)^3/2 = πcosθ^+/8 asinθ^-  dθ^-

using the constraint of equation <ref>

    x dx/(x^2 + {y_o - tanθ^+(x - x_o) }^2)^3/2 = πcosθ^+/8 asinθ^-  dθ^-

defining I ≡ x_otanθ^+ + y_o we can rewrite equation <ref>:

    x dx/(x^2^2θ^+ - 2Itanθ^+ x + I^2)^3/2 = πcosθ^+/8 asinθ^-  dθ^-

which upon integration gives:

    y_o - tanθ^+(x-x_o)/I√(x^2^2θ^+ - 2Itanθ^+ x + I^2) - y_o/I√(x_o^2 + y_o^2) = πcosθ^+/8 a (cosθ^- - cosθ^-_o )

which can be solved for θ^- in terms of x:

    θ = ±cos^-1[ 8 a/Iπcosθ^+{y_o - tanθ^+(x-x_o)/√(x^2^2θ^+ - 2Itanθ^+ x + I^2) - y_o/√(x_o^2 + y_o^2)} + cosθ^-_o]


99
KIM Sangtae Kim , Seppo J. Karrila , Microhydrodynamics: Principles and Selected Applications (book)
Lisa Lisa F. Shatz , Physics of fluids, vol. 16 no. 3, (2004)

