Introduction
Electric fields provide a means of controlling the organization, transport, and mechanical response of particles suspended in fluids. Particle motion is governed by electric forces and torques and by their coupling to the surrounding fluid. In electrorheological suspensions, anisotropic polarization interactions promote chains and columns whose resistance to deformation produces a field-dependent yield stress [3 ] [1, 2]. Their response under flow depends on the competition between electric and hydrodynamic forces and the resulting particle rearrangements [4 , 5 , 6 , 7 , 8 ] ; structural relaxation also contributes to their oscillatory mechanical response [9 ] . Nonuniform fields instead exert a net dielectrophoretic force [10 ] , which has been used to characterize, separate, position, and pattern biological cells and to manipulate smaller biological objects [11 , 12 ] . These applications exploit the dependence of polarization on material properties, field frequency, and the electrical response of the suspending medium.
Electric fields can also drive flow in the surrounding fluid. In induced-charge electro-osmosis, an applied field acts on the diffuse charge that it induces near a polarizable surface, producing fluid motion whose characteristic speed is quadratic in the field at weak forcing [13 , 14 ] . Breaking the symmetry of this flow can produce fluid pumping or propel a suspended particle [15 ] , as demonstrated by the electrophoresis of metallodielectric Janus particles [16 ] . Surface heterogeneity also permits field-controlled assembly and changes in collective organization [17 , 18 , 19 ] . Particle shape and nearby electrodes further modify propulsion, hydrodynamic trapping, and electric orientation or rotation [20 , 21 , 22 ] . Although the charge distributions and forcing protocols differ among these systems, they share a coupling between polarization, interfacial charge transport, and the motion of the surrounding fluid.
Here we focus on leaky dielectric particles, for which finite electrical conduction and charge relaxation are essential. The Taylor–Melcher description treats the bulk phases as electrically neutral, with conduction obeying Ohm’s law, while allowing free charge to accumulate at their interface [23 , 24 ] . The interfacial charge is governed by the balance of conduction currents and surface convection; the associated Maxwell stress couples the electric field to the fluid motion. For deformable drops, tangential electric stresses drive interfacial circulation, which can redistribute surface charge and modify the drop deformation. For rigid particles, electric forces and torques can drive translation and rotation, coupling particle motion to charge relaxation. Rotation redistributes surface charge, while translation through a nonuniform electric field changes the field experienced by the particle. This feedback becomes important when charge relaxation and particle motion occur on comparable time scales [24 , 25 ] .
Quincke rotation is a particularly direct manifestation of this feedback: an initially stationary, isotropic particle spontaneously rotates in a uniform, time-independent electric field [26 , 27 , 28 ] . When charge relaxes more slowly inside the particle than in the fluid, the relaxing part of the induced dipole is antiparallel to the imposed field. A small perturbation tilts the relaxing component of the induced dipole, generating an electric torque that can amplify the perturbation. Above a critical field, the feedback between rotation and surface-charge convection destabilizes the nonrotating state. The resulting steady rotation is sustained by electric torque balancing viscous resistance. In the ideal inertialess problem for an isolated sphere, the rotation axis lies perpendicular to the field but its azimuth is arbitrary. With rotational inertia, the dynamics can be richer: experiments on insulating rotors have demonstrated transitions from steady rotation to chaotic motion, including Lorenz-type chaos [29 , 30 ] .
Quincke rotation also affects suspension rheology and electrical transport. Experiments show that it can reduce the apparent viscosity by assisting an imposed shear [31 ] , while weak shear can reorient the rotation axes and alter the low-shear response [32 ] . The competition between shear, electric attraction, and rotationally generated flow can produce layering and phase separation [33 ] . Viscosity reduction also increases capillary flow rates and changes pressure-driven velocity profiles [34 , 35 ] . Although theory predicts negative apparent viscosity in suitable regimes [36 ] , the pressure-driven experiments found a reduction without reaching negative viscosity [35 ] . Continuum theories coupling particle rotation to the fluid predict polar ordering and spontaneous flows, including flow below the isolated-particle critical field [37 , 38 ] . Quincke rotation also enhances the electrical conductivity of the suspension [39 , 40 ] . These effects connect the rotational instability to suspension microstructure and macroscopic transport in the bulk.
Near a solid boundary, hydrodynamic coupling converts Quincke rotation into translation, producing self-propelled colloidal rollers. These particles continually convert externally supplied electrical energy into motion and have become a model system for investigating active matter. Their motion and interactions can be described through the electric potential and fluid flow, governed by Laplace’s and Stokes equations and coupled through surface-charge dynamics. Populations of spherical rollers form traveling bands and macroscopically directed polar fluids through interactions mediated by the electric field and the surrounding fluid [41 ] . Confinement organizes this motion into persistent vortices [42 ] . Disorder can redirect a flock into sparse channels and eventually destroy polar order [43 ] , whereas obstacles in an annular track can induce repeated flow reversals or arrays of counter-rotating vortices [44 ] .
The electrical driving also provides control over the persistence and collective dynamics of the motion. Pulsing the field allows charge to relax between runs, producing tunable run-and-tumble trajectories and dynamic clustering [45 ] . In larger populations, pulsed driving supports spontaneous shock waves whose fronts propagate faster than the particles they entrain [46 ] . Attraction generated by electrohydrodynamic flow can compete with Quincke propulsion: in experiments on small spherical rollers, initially compact crystallites become motile and continually reshape, merge, and split before dispersing into an active gas as the field increases [47 ] . The field therefore controls several coupled processes, including propulsion, charge memory, and the interactions responsible for aggregation.
Particle shape introduces further ways of converting rotation into transport. Calculations for asymmetric rigid particles show that Quincke rotation can generate propulsion in a bulk fluid, without a nearby wall [48 ] . Experiments with pear-shaped rollers reveal field-dependent curved trajectories and collective states ranging from vortices to rotating flocks [49 ] ; their chiral motion can also suppress long-wavelength density fluctuations while preserving structure at intermediate scales [50 ] . Quincke dumbbells exhibit spinning and orbital motion, and collisions can bind them into rotating tetramers and larger clusters [51 ] . Magnetic control offers an additional independent handle: magnetic Quincke rollers form active chains and can be steered or confined by magnetic fields [52 ] , while other magnetic particles exhibit coexisting rolling directions parallel and perpendicular to the field [53 ] .
Recent experiments also extend the dynamics beyond steady spinning. Small colloids near electrodes can undergo sustained Quincke oscillations, for which ion transport and the spatial variation of conductivity near the electrode are essential [54 ] . Interactions between these oscillators produce simultaneous phase synchronization and alignment of the oscillation axes, with different collective organization in dilute suspensions and dense clusters [55 ] . In a distinct realization with spatially anchored, three-dimensionally printed Quincke micromotors, interactions organize neighboring motors into alternating rotation directions, coherent precession, and propagating phase waves [56 ] . These systems highlight the importance of coupled rotational dynamics, although their oscillation or precession phases are different from the azimuthal orientation of the steady rotating state of an isotropic sphere.
A common ingredient in these diverse collective behaviors is the interaction between electrically driven particles, mediated by both the electric field and the surrounding fluid. Each particle perturbs the local electric field and fluid vorticity experienced by its neighbors, modifying their rotation speeds, selecting their rotation axes, and shifting the critical field for instability. Near a solid boundary, additional electric and hydrodynamic interactions with the surface complicate this coupling. An unbounded fluid therefore provides a simpler starting point for establishing the particle–particle interaction physics.
Previous work developed a far-field model for two Quincke spheres that captured geometry-dependent changes in onset and coupled rotational and translational dynamics [25 ] . Its self-consistent coupled equations provide useful predictions and retain selected higher-order interaction contributions, but these do not constitute a complete asymptotic expansion at those higher orders. We refer to the coupled equations as the unexpanded dipole–rotlet model. Direct numerical simulations have also resolved hydrodynamic interactions between Quincke rollers near a substrate, including lubrication effects [57 ] . Those calculations represent electrical interactions through dipoles and prescribe the Quincke rotation speed from an isolated-particle steady solution.
Here, we study two identical spheres at a fixed separation and orientation in an unbounded fluid subject to a uniform dc electric field. Freezing the pair geometry isolates the local rotational coupling relevant to a suspension. Starting from electric dipole and hydrodynamic rotlet interactions, we develop the first systematic multiple-scale asymptotic theory of the coupled rotational dynamics of two particles undergoing Quincke rotation, covering regimes both near and away from onset. Two effects limit a regular far-field expansion. Because every rotation axis perpendicular to the field is equally possible for an isolated sphere, weak interactions can turn the axes slowly. A correction calculated with fixed axes then grows in time. Independently, the decay rate of rotation-magnitude perturbations vanishes at onset, so the interaction correction can cease to be small even for widely separated particles.
Away from onset, a multiple-scale expansion with a fast and a slow time gives phase equations describing the evolution and selection of the rotation axes. Near onset, we take the distance from the isolated critical field to be of the same order as the interaction strength R − 3 R^{-3} R − 3
We also develop a new boundary-element formulation for two interacting spheres undergoing Quincke rotation. This extends the boundary element method (BEM) previously used for isolated drops [58 ] and isolated rigid particles [48 ] to the coupled two-particle problem. The simulations evolve surface charge on both spheres and resolve their electric and hydrodynamic interactions through the coupled boundary-value problems. They provide a numerical benchmark for the asymptotic theory and the unexpanded dipole–rotlet model at finite separation.
Section II describes the problem set-up, governing equations and dimensionless parameters. Sections III–V develop the regular expansion and slow phase equations away from onset, then establish their near-onset limitation. Sections VI–VIII derive the near-onset amplitude equations and their rotation patterns, critical fields, and interaction mechanisms. Section IX returns to steady speeds away from onset, before Secs. X and XI present the BEM formulation and comparisons.
Problem set-up
We consider two identical rigid spheres of radius a a a E 0 e z E_0\boldsymbol{e}_z E 0 e z Figure 1 . To isolate rotational coupling, we freeze the pair separation and orientation while the surface charge and angular velocities evolve. In a suspension, free translation would also change the pair configuration [59 ] . The hydrodynamic calculation is performed with the electric torques prescribed and the net forces on both spheres set to zero. The resulting translational velocities are discarded. The fluid has permittivity ϵ 1 \epsilon_1 ϵ 1 σ 1 \sigma_1 σ 1 μ \mu μ ϵ 2 \epsilon_2 ϵ 2 σ 2 \sigma_2 σ 2
Figure 1. Geometry of two identical Quincke rotors in a uniform dc electric field E 0 e z E_0\boldsymbol{e}_z E 0 e z x c , 1 \boldsymbol{x}_{c,1} x c , 1 x c , 2 \boldsymbol{x}_{c,2} x c , 2 R R ^ R\hat{\boldsymbol{R}} R R ^ a a a Θ \Theta Θ
The material contrasts and Maxwell–Wagner relaxation time are
ϵ 21 = ϵ 2 − ϵ 1 ϵ 2 + 2 ϵ 1 , σ 21 = σ 2 − σ 1 σ 2 + 2 σ 1 , τ M W = ϵ 2 + 2 ϵ 1 σ 2 + 2 σ 1 . (1) \epsilon_{21} = \frac{\epsilon_2-\epsilon_1}{\epsilon_2+2\epsilon_1}, \quad\sigma_{21} = \frac{\sigma_2-\sigma_1}{\sigma_2+2\sigma_1}, \quad\tau_{\mathrm{MW}} = \frac{\epsilon_2+2\epsilon_1}{\sigma_2+2\sigma_1}. \tag*{(1)} ϵ 21 = ϵ 2 + 2 ϵ 1 ϵ 2 − ϵ 1 , σ 21 = σ 2 + 2 σ 1 σ 2 − σ 1 , τ MW = σ 2 + 2 σ 1 ϵ 2 + 2 ϵ 1 . (1) Lengths, time, electric fields, dipole coefficients, and angular velocities are scaled by a a a τ M W \tau_{\mathrm{MW}} τ MW E 0 E_0 E 0 a 3 E 0 a^3E_0 a 3 E 0 τ M W − 1 \tau_{\mathrm{MW}}^{-1} τ MW − 1 i i i P i \boldsymbol{P}_i P i E i \boldsymbol{E}_i E i Ω i \boldsymbol{\Omega}_i Ω i P i \boldsymbol{P}_i P i P i p h y s = 4 π ϵ 1 a 3 E 0 P i \boldsymbol{P}_i^{\mathrm{phys}}=4\pi\epsilon_1a^3E_0\boldsymbol{P}_i P i phys = 4 π ϵ 1 a 3 E 0 P i
M a = 2 μ τ M W ϵ 1 E 0 2 . (2) \mathrm{Ma}=\frac{2\mu}{\tau_{\mathrm{MW}}\epsilon_1E_0^2}. \tag*{(2)} Ma = τ MW ϵ 1 E 0 2 2 μ . (2) The Mason number compares the viscous stress μ / τ M W \mu/\tau_{\mathrm{MW}} μ / τ MW ϵ 1 E 0 2 \epsilon_1E_0^2 ϵ 1 E 0 2
M a c = ϵ 21 − σ 21 > 0. (3) \mathrm{Ma}_c=\epsilon_{21}-\sigma_{21}>0. \tag*{(3)} Ma c = ϵ 21 − σ 21 > 0. (3) The positivity condition permits Quincke rotation and is equivalent to τ 2 > τ 1 \tau_2>\tau_1 τ 2 > τ 1 τ 1 = ϵ 1 / σ 1 \tau_1=\epsilon_1/\sigma_1 τ 1 = ϵ 1 / σ 1 τ 2 = ϵ 2 / σ 2 \tau_2=\epsilon_2/\sigma_2 τ 2 = ϵ 2 / σ 2 M a < M a c \mathrm{Ma}<\mathrm{Ma}_{c} Ma < Ma c E 0 > E c ( 0 ) E_{0}>E_{c}^{(0)} E 0 > E c ( 0 )
E c ( 0 ) = [ 2 μ ϵ 1 τ M W M a c ] 1 / 2 . (4) E_{c}^{(0)}=\left[\frac{2\mu}{\epsilon_{1}\tau_{\mathrm{MW}}\mathrm{Ma}_{c}}\right]^{1/2}. \tag*{(4)} E c ( 0 ) = [ ϵ 1 τ MW Ma c 2 μ ] 1/2 . (4) This is the critical electric field for Quincke rotation of an isolated sphere [25 ] . Interactions modify the onset field according to the pair separation and orientation, as determined below.
Let x c , i x_{c,i} x c , i i i i
x c , 2 − x c , 1 = R R ^ , ∣ R ^ ∣ = 1. (5) \boldsymbol{x}_{c,2}-\boldsymbol{x}_{c,1}=R\widehat{\boldsymbol{R}},\qquad\lvert\widehat{\boldsymbol{R}}\rvert=1. \tag*{(5)} x c , 2 − x c , 1 = R R , ∣ R ∣ = 1. (5) We consider the far-field regime R ≫ 1 R\gg1 R ≫ 1
δ = R − 3 , Π = I − 3 R ^ R ^ , cos Θ = R ^ ⋅ e z . (6) \delta=R^{-3},\qquad\boldsymbol{\Pi}=\boldsymbol{I}-3\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}},\qquad\cos\Theta=\widehat{\boldsymbol{R}}\cdot\boldsymbol{e}_{z}. \tag*{(6)} δ = R − 3 , Π = I − 3 R R , cos Θ = R ⋅ e z . (6) Here Θ \Theta Θ 0 ≤ Θ ≤ π / 2 0\leq\Theta\leq\pi/2 0 ≤ Θ ≤ π /2
Retaining the electric dipole field and hydrodynamic rotlet disturbance gives
E i = e z − δ Π ⋅ P j , (7) \boldsymbol{E}_{i}=\boldsymbol{e}_{z}-\delta\boldsymbol{\Pi}\cdot\boldsymbol{P}_{j}, \tag*{(7)} E i = e z − δ Π ⋅ P j , (7) Ω i + δ 2 Π ⋅ Ω j = M a − 1 P i × E i , (8) \boldsymbol{\Omega}_{i}+\frac{\delta}{2}\boldsymbol{\Pi}\cdot\boldsymbol{\Omega}_{j}=\mathrm{Ma}^{-1}\boldsymbol{P}_{i}\times\boldsymbol{E}_{i}, \tag*{(8)} Ω i + 2 δ Π ⋅ Ω j = Ma − 1 P i × E i , (8) P ˙ i + ϵ 21 δ Π ⋅ P ˙ j = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) . (9) \dot{\boldsymbol{P}}_{i}+\epsilon_{21}\delta\boldsymbol{\Pi}\cdot\dot{\boldsymbol{P}}_{j}
=\boldsymbol{\Omega}_{i}\times\left(\boldsymbol{P}_{i}-\epsilon_{21}\boldsymbol{E}_{i}\right)-\left(\boldsymbol{P}_{i}-\sigma_{21}\boldsymbol{E}_{i}\right). \tag*{(9)} P ˙ i + ϵ 21 δ Π ⋅ P ˙ j = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) . (9) where i = 1 , 2 i=1,2 i = 1 , 2 j = 3 − i j=3-i j = 3 − i − δ Π ⋅ P j -\delta\boldsymbol{\Pi}\cdot\boldsymbol{P}_{j} − δ Π ⋅ P j j j j i i i − ( δ / 2 ) Π ⋅ Ω j -(\delta/2)\boldsymbol{\Pi}\cdot\boldsymbol{\Omega}_{j} − ( δ /2 ) Π ⋅ Ω j Π \boldsymbol{\Pi} Π − 2 -2 − 2 1 1 1 [60 ] ; the dipole–rotlet description follows the far-field framework of Ref. [25 ] .
We call Eqs. (7) the unexpanded dipole–rotlet model . Solving the coupled equations without expanding in δ \delta δ δ \delta δ
Differentiating Eq. (7a) gives E ˙ i = − δ Π ⋅ P ˙ j \dot{\boldsymbol{E}}_{i}=-\delta\boldsymbol{\Pi}\cdot\dot{\boldsymbol{P}}_{j} E ˙ i = − δ Π ⋅ P ˙ j ϵ 21 E ˙ i \epsilon_{21}\dot{\boldsymbol{E}}_{i} ϵ 21 E ˙ i
Regular expansion in the interaction strength
At fixed M a \mathrm{Ma} Ma
P i = P i ( 0 ) + δ P i ( 1 ) + O ( δ 2 ) , (10) \boldsymbol{P}_{i}=\boldsymbol{P}_{i}^{(0)}+\delta\boldsymbol{P}_{i}^{(1)}+\mathcal{O}(\delta^{2}), \tag*{(10)} P i = P i ( 0 ) + δ P i ( 1 ) + O ( δ 2 ) , (10) with analogous expansions for E i \boldsymbol{E}_{i} E i Ω i \boldsymbol{\Omega}_{i} Ω i § 1.4 . We substitute these expansions into Eqs. (7) and equate coefficients of equal powers of δ \delta δ
E i ( 0 ) = e z , E i ( 1 ) = − Π ⋅ P j ( 0 ) . (11) \boldsymbol{E}_{i}^{(0)}=\boldsymbol{e}_{z},\qquad\boldsymbol{E}_{i}^{(1)}=-\boldsymbol{\Pi}\cdot\boldsymbol{P}_{j}^{(0)}. \tag*{(11)} E i ( 0 ) = e z , E i ( 1 ) = − Π ⋅ P j ( 0 ) . (11) Ω i ( 0 ) = M a − 1 P i ( 0 ) × e z . (12) \boldsymbol{\Omega}_{i}^{(0)}=\mathrm{Ma}^{-1}\boldsymbol{P}_{i}^{(0)}\times\boldsymbol{e}_{z}. \tag*{(12)} Ω i ( 0 ) = Ma − 1 P i ( 0 ) × e z . (12) Ω i ( 1 ) = M a − 1 [ P i ( 1 ) × e z − P i ( 0 ) × ( Π ⋅ P j ( 0 ) ) − 1 2 Π ⋅ ( P j ( 0 ) × e z ) ] . (13) \boldsymbol{\Omega}_{i}^{(1)}
=\mathrm{Ma}^{-1}\left[\boldsymbol{P}_{i}^{(1)}\times\boldsymbol{e}_{z}
-\boldsymbol{P}_{i}^{(0)}\times\left(\boldsymbol{\Pi}\cdot\boldsymbol{P}_{j}^{(0)}\right)
-\frac{1}{2}\boldsymbol{\Pi}\cdot\left(\boldsymbol{P}_{j}^{(0)}\times\boldsymbol{e}_{z}\right)\right]. \tag*{(13)} Ω i ( 1 ) = Ma − 1 [ P i ( 1 ) × e z − P i ( 0 ) × ( Π ⋅ P j ( 0 ) ) − 2 1 Π ⋅ ( P j ( 0 ) × e z ) ] . (13) The leading dipoles obey the isolated-sphere equations
P ˙ i ( 0 ) = Ω i ( 0 ) × ( P i ( 0 ) − ϵ 21 e z ) − ( P i ( 0 ) − σ 21 e z ) . (14) \dot{\boldsymbol{P}}_{i}^{(0)}
=\boldsymbol{\Omega}_{i}^{(0)}\times\left(\boldsymbol{P}_{i}^{(0)}-\epsilon_{21}\boldsymbol{e}_{z}\right)
-\left(\boldsymbol{P}_{i}^{(0)}-\sigma_{21}\boldsymbol{e}_{z}\right). \tag*{(14)} P ˙ i ( 0 ) = Ω i ( 0 ) × ( P i ( 0 ) − ϵ 21 e z ) − ( P i ( 0 ) − σ 21 e z ) . (14) At order δ \delta δ
P ˙ i ( 1 ) = Ω i ( 1 ) × ( P i ( 0 ) − ϵ 21 e z ) + Ω i ( 0 ) × [ P i ( 1 ) + ϵ 21 Π ⋅ P j ( 0 ) ] − P i ( 1 ) − σ 21 Π ⋅ P j ( 0 ) − ϵ 21 Π ⋅ P ˙ j ( 0 ) . (15) \dot{\boldsymbol{P}}_{i}^{(1)}
=\boldsymbol{\Omega}_{i}^{(1)}\times\left(\boldsymbol{P}_{i}^{(0)}-\epsilon_{21}\boldsymbol{e}_{z}\right)
+\boldsymbol{\Omega}_{i}^{(0)}\times\left[\boldsymbol{P}_{i}^{(1)}+\epsilon_{21}\boldsymbol{\Pi}\cdot\boldsymbol{P}_{j}^{(0)}\right]
-\boldsymbol{P}_{i}^{(1)}-\sigma_{21}\boldsymbol{\Pi}\cdot\boldsymbol{P}_{j}^{(0)}-\epsilon_{21}\boldsymbol{\Pi}\cdot\dot{\boldsymbol{P}}_{j}^{(0)}. \tag*{(15)} P ˙ i ( 1 ) = Ω i ( 1 ) × ( P i ( 0 ) − ϵ 21 e z ) + Ω i ( 0 ) × [ P i ( 1 ) + ϵ 21 Π ⋅ P j ( 0 ) ] − P i ( 1 ) − σ 21 Π ⋅ P j ( 0 ) − ϵ 21 Π ⋅ P ˙ j ( 0 ) . (15) The equations are solved in succession: first Eq. (14) for P i ( 0 ) \boldsymbol{P}_{i}^{(0)} P i ( 0 ) (15) for P i ( 1 ) \boldsymbol{P}_{i}^{(1)} P i ( 1 ) P i ( 0 ) + δ P i ( 1 ) \boldsymbol{P}_{i}^{(0)}+\delta\boldsymbol{P}_{i}^{(1)} P i ( 0 ) + δ P i ( 1 ) δ P i ( 1 ) \delta\boldsymbol{P}_{i}^{(1)} δ P i ( 1 )
Steady rotation and slow changes of the rotation axes
Non-rotating state and isolated rotating family
The labels ‘nr’ and ‘rot’ denote non-rotating and rotating equilibria, respectively; superscripts (0) and (1) continue to denote regular expansion coefficients. The isolated non-rotating state is P n r ( 0 ) = σ 21 e z P_{\mathrm{nr}}^{(0)}=\sigma_{21}\boldsymbol{e}_z P nr ( 0 ) = σ 21 e z Ω i = 0 \boldsymbol{\Omega}_i=\boldsymbol{0} Ω i = 0 P ˙ i = 0 \dot{\boldsymbol{P}}_i=\boldsymbol{0} P ˙ i = 0 P i = σ 21 E i \boldsymbol{P}_i=\sigma_{21}\boldsymbol{E}_i P i = σ 21 E i δ \delta δ
P n r = σ 21 e z − δ σ 21 2 Π ⋅ e z + O ( δ 2 ) , (16) \boldsymbol{P}_{\mathrm{nr}}=\sigma_{21}\boldsymbol{e}_z-\delta\sigma_{21}^{2}\boldsymbol{\Pi}\cdot\boldsymbol{e}_z+O(\delta^2), \tag*{(16)} P nr = σ 21 e z − δ σ 21 2 Π ⋅ e z + O ( δ 2 ) , (16) P n r ( 0 ) = σ 21 e z , P n r ( 1 ) = − σ 21 2 Π ⋅ e z . (17) \boldsymbol{P}_{\mathrm{nr}}^{(0)}=\sigma_{21}\boldsymbol{e}_z,\qquad\boldsymbol{P}_{\mathrm{nr}}^{(1)}=-\sigma_{21}^{2}\boldsymbol{\Pi}\cdot\boldsymbol{e}_z. \tag*{(17)} P nr ( 0 ) = σ 21 e z , P nr ( 1 ) = − σ 21 2 Π ⋅ e z . (17) In the non-rotating state, the dipole is collinear with its local field, so the electric torque vanishes. Their common direction is generally tilted relative to the imposed field. Throughout the paper, transverse means perpendicular to the imposed field E 0 e z E_0\boldsymbol{e}_z E 0 e z P i , ⊥ = ( I − e z e z ) ⋅ P i \boldsymbol{P}_{i,\perp}=(\boldsymbol{I}-\boldsymbol{e}_z\boldsymbol{e}_z)\cdot\boldsymbol{P}_i P i , ⊥ = ( I − e z e z ) ⋅ P i P ( 0 ) = P z ( 0 ) e z + P ⊥ ( 0 ) \boldsymbol{P}^{(0)}=P_z^{(0)}\boldsymbol{e}_z+\boldsymbol{P}_{\perp}^{(0)} P ( 0 ) = P z ( 0 ) e z + P ⊥ ( 0 ) (14) becomes
P ˙ ⊥ ( 0 ) = [ λ 0 − P z ( 0 ) − σ 21 M a ] P ⊥ ( 0 ) , (18) \dot{\boldsymbol{P}}_{\perp}^{(0)}=\left[\lambda_0-\frac{P_z^{(0)}-\sigma_{21}}{\mathrm{Ma}}\right]\boldsymbol{P}_{\perp}^{(0)}, \tag*{(18)} P ˙ ⊥ ( 0 ) = [ λ 0 − Ma P z ( 0 ) − σ 21 ] P ⊥ ( 0 ) , (18) P ˙ z ( 0 ) = M a − 1 ∣ P ⊥ ( 0 ) ∣ 2 − ( P z ( 0 ) − σ 21 ) , (19) \dot{P}_z^{(0)}=\mathrm{Ma}^{-1}\lvert\boldsymbol{P}_{\perp}^{(0)}\rvert^2-\left(P_z^{(0)}-\sigma_{21}\right), \tag*{(19)} P ˙ z ( 0 ) = Ma − 1 ∣ P ⊥ ( 0 ) ∣ 2 − ( P z ( 0 ) − σ 21 ) , (19) where
λ 0 = M a c M a − 1 = ( E 0 E c ( 0 ) ) 2 − 1. (20) \lambda_0=\frac{\mathrm{Ma}_c}{\mathrm{Ma}}-1=\left(\frac{E_0}{E_c^{(0)}}\right)^2-1. \tag*{(20)} λ 0 = Ma Ma c − 1 = ( E c ( 0 ) E 0 ) 2 − 1. (20) A small transverse dipole perturbation about the non-rotating state grows or decays at rate λ 0 \lambda_0 λ 0 λ 0 > 0 \lambda_0>0 λ 0 > 0
P r o t ( 0 ) ( ϕ ) = ( ϵ 21 − M a ) e z + P ⊥ ( 0 ) e ρ ( ϕ ) , \boldsymbol{P}_{\mathrm{rot}}^{(0)}(\phi)=(\epsilon_{21}-\mathrm{Ma})\boldsymbol{e}_z+P_{\perp}^{(0)}\boldsymbol{e}_{\rho}(\phi), P rot ( 0 ) ( ϕ ) = ( ϵ 21 − Ma ) e z + P ⊥ ( 0 ) e ρ ( ϕ ) , P ⊥ ( 0 ) = M a λ 0 . (21) P_{\perp}^{(0)}=\mathrm{Ma}\sqrt{\lambda_0}. \tag*{(21)} P ⊥ ( 0 ) = Ma λ 0 . (21) Here P ⊥ ( 0 ) P_{\perp}^{(0)} P ⊥ ( 0 ) ∣ P ⊥ ( 0 ) ∣ \lvert\boldsymbol{P}_{\perp}^{(0)}\rvert ∣ P ⊥ ( 0 ) ∣
e ρ = ( cos ϕ , sin ϕ , 0 ) , e ϕ = e z × e ρ . (22) \boldsymbol{e}_{\rho}=(\cos\phi,\sin\phi,0),\qquad\boldsymbol{e}_{\phi}=\boldsymbol{e}_z\times\boldsymbol{e}_{\rho}. \tag*{(22)} e ρ = ( cos ϕ , sin ϕ , 0 ) , e ϕ = e z × e ρ . (22) These are the radial and azimuthal unit vectors in the plane perpendicular to the imposed field, with e ρ \boldsymbol{e}_{\rho} e ρ ϕ \phi ϕ (21) into Eq. (12) gives Ω r o t ( 0 ) = − λ 0 e ϕ \boldsymbol{\Omega}_{\mathrm{rot}}^{(0)}=-\sqrt{\lambda_0}\boldsymbol{e}_{\phi} Ω rot ( 0 ) = − λ 0 e ϕ λ 0 \sqrt{\lambda_0} λ 0
Solvability of the first-order steady problem
We now choose the leading dipoles to be steady rotating solutions (21) , with prescribed azimuthal angles ϕ i \phi_i ϕ i
We use phase to mean this azimuthal angle throughout the analysis. Since P ˙ j ( 0 ) = 0 \dot{\boldsymbol{P}}_j^{(0)}=\boldsymbol{0} P ˙ j ( 0 ) = 0 (15) vanishes for this choice of leading solution. Substituting Eq. (13) into Eq. (15) and collecting the terms proportional to P i ( 1 ) \boldsymbol{P}_i^{(1)} P i ( 1 )
P ˙ i ( 1 ) = J i ⋅ P i ( 1 ) + G i . (23) \dot{\boldsymbol{P}}_i^{(1)}=\boldsymbol{J}_i\cdot\boldsymbol{P}_i^{(1)}+\boldsymbol{G}_i. \tag*{(23)} P ˙ i ( 1 ) = J i ⋅ P i ( 1 ) + G i . (23) Here J i \boldsymbol{J}_i J i P i ( 0 ) \boldsymbol{P}_i^{(0)} P i ( 0 ) G i \boldsymbol{G}_i G i Ω i ( 1 ) \boldsymbol{\Omega}_i^{(1)} Ω i ( 1 ) (13) produced by particle j j j
Ω i , i n t ( 1 ) = − M a − 1 [ P i ( 0 ) × ( Π ⋅ P j ( 0 ) ) + 1 2 Π ⋅ ( P j ( 0 ) × e z ) ] . (24) \boldsymbol{\Omega}_{i,\mathrm{int}}^{(1)}=-\mathrm{Ma}^{-1}\left[\boldsymbol{P}_i^{(0)}\times(\boldsymbol{\Pi}\cdot\boldsymbol{P}_j^{(0)})+\frac{1}{2}\boldsymbol{\Pi}\cdot(\boldsymbol{P}_j^{(0)}\times\boldsymbol{e}_z)\right]. \tag*{(24)} Ω i , int ( 1 ) = − Ma − 1 [ P i ( 0 ) × ( Π ⋅ P j ( 0 ) ) + 2 1 Π ⋅ ( P j ( 0 ) × e z ) ] . (24) Substitution into Eq. (15) gives
G i = Ω i , i n t ( 1 ) × ( P i ( 0 ) − ϵ 21 e z ) + ϵ 21 Ω i ( 0 ) × ( Π ⋅ P j ( 0 ) ) − σ 21 Π ⋅ P j ( 0 ) . (25) \boldsymbol{G}_i=\boldsymbol{\Omega}_{i,\mathrm{int}}^{(1)}\times(\boldsymbol{P}_i^{(0)}-\epsilon_{21}\boldsymbol{e}_z)+\epsilon_{21}\boldsymbol{\Omega}_i^{(0)}\times(\boldsymbol{\Pi}\cdot\boldsymbol{P}_j^{(0)})-\sigma_{21}\boldsymbol{\Pi}\cdot\boldsymbol{P}_j^{(0)}. \tag*{(25)} G i = Ω i , int ( 1 ) × ( P i ( 0 ) − ϵ 21 e z ) + ϵ 21 Ω i ( 0 ) × ( Π ⋅ P j ( 0 ) ) − σ 21 Π ⋅ P j ( 0 ) . (25) Resolve the correction along e ρ i \boldsymbol{e}_{\rho_i} e ρ i e ϕ i \boldsymbol{e}_{\phi_i} e ϕ i e z \boldsymbol{e}_z e z
[ J i ] ( ρ , ϕ , z ) = ( 0 0 − P ⊥ ( 0 ) / M a 0 0 0 2 P ⊥ ( 0 ) / M a 0 − 1 ) . (26) [\boldsymbol{J}_i]_{(\rho,\phi,z)}=
\begin{pmatrix}
0 & 0 & -P_{\perp}^{(0)}/\mathrm{Ma}\\
0 & 0 & 0\\
2P_{\perp}^{(0)}/\mathrm{Ma} & 0 & -1
\end{pmatrix}. \tag*{(26)} [ J i ] ( ρ , ϕ , z ) = 0 0 2 P ⊥ ( 0 ) / Ma 0 0 0 − P ⊥ ( 0 ) / Ma 0 − 1 . (26) The second row is zero: a small change of the isolated sphere’s phase produces another steady state, so there is no restoring term for that change. The three component equations are therefore
P ˙ i , ρ ( 1 ) = − P ⊥ ( 0 ) M a P i , z ( 1 ) + G i , ρ , \dot{P}_{i,\rho}^{(1)}=-\frac{P_{\perp}^{(0)}}{\mathrm{Ma}}P_{i,z}^{(1)}+G_{i,\rho}, P ˙ i , ρ ( 1 ) = − Ma P ⊥ ( 0 ) P i , z ( 1 ) + G i , ρ , P ˙ i , ϕ ( 1 ) = G i , ϕ , \dot{P}_{i,\phi}^{(1)}=G_{i,\phi}, P ˙ i , ϕ ( 1 ) = G i , ϕ , P ˙ i , z ( 1 ) = 2 P ⊥ ( 0 ) M a P i , ρ ( 1 ) − P i , z ( 1 ) + G i , z . (27) \dot{P}_{i,z}^{(1)}=\frac{2P_{\perp}^{(0)}}{\mathrm{Ma}}P_{i,\rho}^{(1)}-P_{i,z}^{(1)}+G_{i,z}. \tag*{(27)} P ˙ i , z ( 1 ) = Ma 2 P ⊥ ( 0 ) P i , ρ ( 1 ) − P i , z ( 1 ) + G i , z . (27) For a time-independent first-order correction to steady rotation, all three left-hand sides must vanish. The second equation has a steady solution only if its forcing vanishes. Evaluating the e ϕ i \boldsymbol{e}_{\phi_i} e ϕ i (25) gives
G i , ϕ = P ⊥ ( 0 ) 2 ( e ρ i + P ⊥ ( 0 ) M a e z ) ⋅ Π ⋅ e ϕ j = 0 , i = 1 , 2. (28) \begin{aligned}
G_{i,\phi}=\frac{P_{\perp}^{(0)}}{2}\left(\boldsymbol{e}_{\rho_i}+\frac{P_{\perp}^{(0)}}{\mathrm{Ma}}\boldsymbol{e}_z\right)\cdot\boldsymbol{\Pi}\cdot\boldsymbol{e}_{\phi_j} \\
&=0,\qquad i=1,2. \tag*{(28)}
\end{aligned} G i , ϕ = 2 P ⊥ ( 0 ) ( e ρ i + Ma P ⊥ ( 0 ) e z ) ⋅ Π ⋅ e ϕ j = 0 , i = 1 , 2. (28) This is the solvability condition for the first-order steady problem. Equivalently, multiplying J i P i ( 1 ) = − G i \boldsymbol{J}_i\boldsymbol{P}_i^{(1)}=-\boldsymbol{G}_i J i P i ( 1 ) = − G i e ϕ i \boldsymbol{e}_{\phi_i} e ϕ i [61 ] . When this condition holds, the first and third equations give
P i , z ( 1 ) = M a P ⊥ ( 0 ) G i , ρ , (29) P_{i,z}^{(1)} = \frac{\mathrm{Ma}}{P_{\perp}^{(0)}}G_{i,\rho}, \tag*{(29)} P i , z ( 1 ) = P ⊥ ( 0 ) Ma G i , ρ , (29) P i , ρ ( 1 ) = M a 2 P ⊥ ( 0 ) ( M a P ⊥ ( 0 ) G i , ρ − G i , z ) . (30) P_{i,\rho}^{(1)} = \frac{\mathrm{Ma}}{2P_{\perp}^{(0)}}\left(\frac{\mathrm{Ma}}{P_{\perp}^{(0)}}G_{i,\rho}-G_{i,z}\right). \tag*{(30)} P i , ρ ( 1 ) = 2 P ⊥ ( 0 ) Ma ( P ⊥ ( 0 ) Ma G i , ρ − G i , z ) . (30) A constant P i , ϕ ( 1 ) P_{i,\phi}^{(1)} P i , ϕ ( 1 ) δ \delta δ
Secular growth and a slow time for the phases
For fixed azimuthal angles ϕ 1 \phi_1 ϕ 1 ϕ 2 \phi_2 ϕ 2 (28) , the azimuthal component of Eq. (27) gives
P i , ϕ ( 1 ) ( t ) = P i , ϕ ( 1 ) ( 0 ) + G i , ϕ t . (31) P_{i,\phi}^{(1)}(t)=P_{i,\phi}^{(1)}(0)+G_{i,\phi}t. \tag*{(31)} P i , ϕ ( 1 ) ( t ) = P i , ϕ ( 1 ) ( 0 ) + G i , ϕ t . (31) The term proportional to t t t secular term : it eventually makes the dipole correction δ P i , ϕ ( 1 ) \delta P_{i,\phi}^{(1)} δ P i , ϕ ( 1 ) P ⊥ ( 0 ) P_{\perp}^{(0)} P ⊥ ( 0 )
e ρ ( ϕ i + Δ ϕ i ) = e ρ ( ϕ i ) + Δ ϕ i e ϕ i + O ( Δ ϕ i 2 ) , (32) \boldsymbol{e}_{\rho}(\phi_i+\Delta\phi_i)=\boldsymbol{e}_{\rho}(\phi_i)+\Delta\phi_i\boldsymbol{e}_{\phi_i}+O(\Delta\phi_i^2), \tag*{(32)} e ρ ( ϕ i + Δ ϕ i ) = e ρ ( ϕ i ) + Δ ϕ i e ϕ i + O ( Δ ϕ i 2 ) , (32) so
Δ ϕ i ≃ δ P i , ϕ ( 1 ) P ⊥ ( 0 ) . (33) \Delta\phi_i \simeq\frac{\delta P_{i,\phi}^{(1)}}{P_{\perp}^{(0)}}. \tag*{(33)} Δ ϕ i ≃ P ⊥ ( 0 ) δ P i , ϕ ( 1 ) . (33) To remove the secular growth, we continue to expand about the isolated rotating family in Eq. (21) , with λ 0 > 0 \lambda_0>0 λ 0 > 0 δ → 0 \delta\to0 δ → 0 P ⊥ ( 0 ) = M a λ 0 P_{\perp}^{(0)}=\mathrm{Ma}\sqrt{\lambda_0} P ⊥ ( 0 ) = Ma λ 0 [61 , 62 ] , retaining the original time t t t τ = δ t \tau=\delta t τ = δ t τ \tau τ
P i = P r o t ( 0 ) ( ϕ i ( τ ) ) + δ P i ( 1 ) ( t , τ ) + ⋯ . (34) \boldsymbol{P}_i=\boldsymbol{P}_{\mathrm{rot}}^{(0)}(\phi_i(\tau))+\delta\boldsymbol{P}_i^{(1)}(t,\tau)+\cdots. \tag*{(34)} P i = P rot ( 0 ) ( ϕ i ( τ )) + δ P i ( 1 ) ( t , τ ) + ⋯ . (34) Here P i ( 1 ) ( t , τ ) \boldsymbol{P}_i^{(1)}(t,\tau) P i ( 1 ) ( t , τ ) t t t
During the calculation, t t t τ \tau τ
d d t = ∂ t + δ ∂ τ . (35) \frac{\mathrm{d}}{\mathrm{d}t}=\partial_t+\delta\partial_\tau. \tag*{(35)} d t d = ∂ t + δ ∂ τ . (35) Substitution into Eqs. (7) gives, at order δ \delta δ
∂ t P i ( 1 ) = J i ⋅ P i ( 1 ) + G i − P i ( 0 ) d ϕ i d τ e ϕ i . (36) \partial_t\boldsymbol{P}_i^{(1)}=\boldsymbol{J}_i\cdot\boldsymbol{P}_i^{(1)}+\boldsymbol{G}_i-P_i^{(0)}\frac{\mathrm{d}\phi_i}{\mathrm{d}\tau}\boldsymbol{e}_{\phi_i}. \tag*{(36)} ∂ t P i ( 1 ) = J i ⋅ P i ( 1 ) + G i − P i ( 0 ) d τ d ϕ i e ϕ i . (36) Its component along e ϕ i \boldsymbol{e}_{\phi_i} e ϕ i
P ( 0 ) d ϕ i d τ = G i , ϕ ( ϕ 1 , ϕ 2 ) . (37) P^{(0)}\frac{\mathrm{d}\phi_i}{\mathrm{d}\tau}=G_{i,\phi}(\phi_1,\phi_2). \tag*{(37)} P ( 0 ) d τ d ϕ i = G i , ϕ ( ϕ 1 , ϕ 2 ) . (37) Equation (37) determines the slow evolution of the rotation axes; its equilibria satisfy G i , ϕ = 0 G_{i,\phi}=0 G i , ϕ = 0
Phase selection away from onset
We now determine which rotating states are stable to small changes in the directions of their rotation axes. We consider λ 0 > 0 \lambda_0>0 λ 0 > 0 δ / λ 0 ≪ 1 \delta/\lambda_0\ll1 δ / λ 0 ≪ 1
Choose the line of centres in the x x x z z z R ^ = ( sin Θ , 0 , cos Θ ) \hat{\boldsymbol{R}}=(\sin\Theta,0,\cos\Theta) R ^ = ( sin Θ , 0 , cos Θ ) (37) then becomes
d ϕ i d τ = 1 2 [ sin ϕ i cos ϕ j − ( Π x x cos ϕ i + λ 0 Π z x ) sin ϕ j ] . (38) \frac{\mathrm{d}\phi_i}{\mathrm{d}\tau}=\frac{1}{2}\left[\sin\phi_i\cos\phi_j-\left(\Pi_{xx}\cos\phi_i+\sqrt{\lambda_0}\Pi_{zx}\right)\sin\phi_j\right]. \tag*{(38)} d τ d ϕ i = 2 1 [ sin ϕ i cos ϕ j − ( Π xx cos ϕ i + λ 0 Π z x ) sin ϕ j ] . (38) where Π x x = 3 cos 2 Θ − 2 \Pi_{xx}=3\cos^2\Theta-2 Π xx = 3 cos 2 Θ − 2 Π z x = − 3 sin Θ cos Θ \Pi_{zx}=-3\sin\Theta\cos\Theta Π z x = − 3 sin Θ cos Θ (21) and (12) , Ω i = − λ 0 e ϕ i \boldsymbol{\Omega}_i=-\sqrt{\lambda_0}\,\boldsymbol{e}_{\phi_i} Ω i = − λ 0 e ϕ i π \pi π
We first examine the counter-rotating branch ( ϕ 1 , ϕ 2 ) = ( 0 , π ) ; ( π , 0 ) (\phi_1,\phi_2)=(0,\pi);(\pi,0) ( ϕ 1 , ϕ 2 ) = ( 0 , π ) ; ( π , 0 ) x x x y y y (38) makes both evolution rates vanish for every pair angle. This branch therefore exists throughout the parameter range considered.
To determine its stability, we linearize Eq. (38) with respect to the two phases. The resulting phase Jacobian has trace − 1 -1 − 1
det J ϕ = 3 4 sin 2 Θ [ 3 ( 1 + λ 0 ) cos 2 Θ − 1 ] . (39) \det J_{\phi}=\frac{3}{4}\sin^2\Theta\left[3(1+\lambda_0)\cos^2\Theta-1\right]. \tag*{(39)} det J ϕ = 4 3 sin 2 Θ [ 3 ( 1 + λ 0 ) cos 2 Θ − 1 ] . (39) For 0 < Θ < π / 2 0<\Theta<\pi/2 0 < Θ < π /2 Θ < Θ ∗ \Theta<\Theta_* Θ < Θ ∗ Θ > Θ ∗ \Theta> \Theta_\ast Θ > Θ ∗
cos 2 Θ ∗ = 1 3 ( 1 + λ 0 ) , Θ ∗ = cos − 1 ( E c ( 0 ) 3 E 0 ) . (40) \cos^2 \Theta_\ast= \frac{1}{3(1+\lambda_0)}, \qquad\Theta_\ast= \cos^{-1}\left(\frac{E_{\mathrm{c}}^{(0)}}{\sqrt{3}E_0}\right). \tag*{(40)} cos 2 Θ ∗ = 3 ( 1 + λ 0 ) 1 , Θ ∗ = cos − 1 ( 3 E 0 E c ( 0 ) ) . (40) At Θ = 0 \Theta= 0 Θ = 0
We next examine co-rotating states. Setting ϕ 1 = ϕ 2 = ϕ \phi_1 = \phi_2 = \phi ϕ 1 = ϕ 2 = ϕ (38) , the equilibrium condition for co-rotating states factorizes as
0 = 1 2 sin ϕ [ ( 1 − Π x x ) cos ϕ − λ 0 Π x z ] . (41) 0 = \frac{1}{2}\sin\phi\left[(1-\Pi_{xx})\cos\phi-\sqrt{\lambda_0}\Pi_{xz}\right]. \tag*{(41)} 0 = 2 1 sin ϕ [ ( 1 − Π xx ) cos ϕ − λ 0 Π x z ] . (41) There are consequently two types of co-rotating solution. The first has ϕ = 0 \phi= 0 ϕ = 0 π \pi π x x x + 1 +1 + 1
The second co-rotating family satisfies
ϕ 1 = ϕ 2 = ϕ , cos ϕ = − λ 0 cot Θ . (42) \phi_1 = \phi_2 = \phi, \qquad\cos\phi= -\sqrt{\lambda_0}\cot\Theta. \tag*{(42)} ϕ 1 = ϕ 2 = ϕ , cos ϕ = − λ 0 cot Θ. (42) Here the common dipole direction changes with pair angle and field strength. A real phase requires ∣ cos ϕ ∣ ≤ 1 |\cos\phi| \le1 ∣ cos ϕ ∣ ≤ 1
Θ ≥ Θ e x , Θ e x = tan − 1 λ 0 = cos − 1 ( E c ( 0 ) E 0 ) . (43) \Theta\ge\Theta_{\mathrm{ex}}, \qquad\Theta_{\mathrm{ex}} = \tan^{-1}\sqrt{\lambda_0} = \cos^{-1}\left(\frac{E_{\mathrm{c}}^{(0)}}{E_0}\right). \tag*{(43)} Θ ≥ Θ ex , Θ ex = tan − 1 λ 0 = cos − 1 ( E 0 E c ( 0 ) ) . (43) For Θ > Θ e x \Theta> \Theta_{\mathrm{ex}} Θ > Θ ex sin ϕ ≠ 0 \sin\phi\ne0 sin ϕ = 0 Θ = Θ e x \Theta= \Theta_{\mathrm{ex}} Θ = Θ ex cos ϕ = − 1 \cos\phi= -1 cos ϕ = − 1 ( ϕ 1 , ϕ 2 ) = ( π , π ) (\phi_1,\phi_2) = (\pi,\pi) ( ϕ 1 , ϕ 2 ) = ( π , π ) (42) , perturbations that change both phases equally have growth rate ν s \nu_s ν s ν a \nu_a ν a
ν s = − 3 2 sin 2 Θ sin 2 ϕ , (44) \nu_s = -\frac{3}{2}\sin^2\Theta\sin^2\phi, \tag*{(44)} ν s = − 2 3 sin 2 Θ sin 2 ϕ , (44) ν a = 1 2 [ 3 ( 1 + λ 0 ) cos 2 Θ − 1 ] . (45) \nu_a = \frac{1}{2}\left[3(1+\lambda_0)\cos^2\Theta-1\right]. \tag*{(45)} ν a = 2 1 [ 3 ( 1 + λ 0 ) cos 2 Θ − 1 ] . (45) These rates are measured on the slow time τ \tau τ ν s < 0 \nu_s < 0 ν s < 0 ν a \nu_a ν a Θ < Θ ∗ \Theta< \Theta_\ast Θ < Θ ∗ Θ > Θ ∗ \Theta> \Theta_\ast Θ > Θ ∗ Θ e x < Θ < Θ ∗ \Theta_{\mathrm{ex}} < \Theta< \Theta_\ast Θ ex < Θ < Θ ∗ Θ > Θ ∗ \Theta> \Theta_\ast Θ > Θ ∗ Θ = π / 2 \Theta= \pi/2 Θ = π /2 e y \mathbf e_y e y − e y -\mathbf e_y − e y x x x ϕ 2 = − ϕ 1 \phi_2 = -\phi_1 ϕ 2 = − ϕ 1
As an illustration, consider a line of centres perpendicular to the imposed field, Θ = π / 2 \Theta= \pi/2 Θ = π /2 R ^ = e x \hat{\mathbf R} = \mathbf e_x R ^ = e x i i i
− δ 2 Π ⋅ Ω j ( 0 ) = δ Ω j , x ( 0 ) e x − δ 2 Ω j , y ( 0 ) e y . (46) -\frac{\delta}{2}\boldsymbol{\Pi}\cdot\boldsymbol{\Omega}_j^{(0)}
= \delta\Omega_{j,x}^{(0)}\mathbf e_x-\frac{\delta}{2}\Omega_{j,y}^{(0)}\mathbf e_y. \tag*{(46)} − 2 δ Π ⋅ Ω j ( 0 ) = δ Ω j , x ( 0 ) e x − 2 δ Ω j , y ( 0 ) e y . (46) where Ω j , z ( 0 ) = 0 \Omega_{j,z}^{(0)}=0 Ω j , z ( 0 ) = 0 y y y ( 0 , 0 ) (0,0) ( 0 , 0 ) ( π , π ) (\pi,\pi) ( π , π ) y y y x x x ϕ = π / 2 \phi=\pi/2 ϕ = π /2 3 π / 2 3\pi/2 3 π /2 y y y
For perturbations that tilt both axes together, the equal-phase equation at Θ = π / 2 \Theta=\pi/2 Θ = π /2
d ϕ d τ = 3 4 sin ( 2 ϕ ) . (47) \frac{d\phi}{d\tau}=\frac{3}{4}\sin(2\phi). \tag*{(47)} d τ d ϕ = 4 3 sin ( 2 ϕ ) . (47) Its linear growth rate is + 3 / 2 +3/2 + 3/2 ϕ = 0 \phi=0 ϕ = 0 π \pi π − 3 / 2 -3/2 − 3/2 ϕ = π / 2 \phi=\pi/2 ϕ = π /2 3 π / 2 3\pi/2 3 π /2 ν a = − 1 / 2 \nu_a=-1/2 ν a = − 1/2
Figure 2 distinguishes the existence boundary Θ e x \Theta_{\mathrm{ex}} Θ ex Θ ∗ \Theta_\ast Θ ∗ E 0 / E c ( 0 ) = 2 E_0/E_{\mathrm{c}}^{(0)}=2 E 0 / E c ( 0 ) = 2 λ 0 = 3 \lambda_0=3 λ 0 = 3 Θ e x = 60 ∘ \Theta_{\mathrm{ex}}=60^\circ Θ ex = 6 0 ∘ Θ ∗ = 73.2 ∘ \Theta_\ast=73.2^\circ Θ ∗ = 73. 2 ∘ (42) exist, but only counter-rotation is stable. Above Θ ∗ \Theta_* Θ ∗ ϕ = 0 \phi=0 ϕ = 0 π \pi π
Figure 2. Phase selection away from onset. The solid stability boundary Θ ∗ \Theta_\ast Θ ∗ ( ϕ 1 , ϕ 2 ) = ( 0 , π ) (\phi_1,\phi_2)=(0,\pi) ( ϕ 1 , ϕ 2 ) = ( 0 , π ) (42) above it. That family exists above the dashed boundary Θ e x \Theta_{\mathrm{ex}} Θ ex ( 0 , 0 ) (0,0) ( 0 , 0 ) ( π , π ) (\pi,\pi) ( π , π ) Θ ∗ = 73.2 ∘ \Theta_\ast=73.2^\circ Θ ∗ = 73. 2 ∘ E 0 / E c ( 0 ) = 2 E_0/E_{\mathrm{c}}^{(0)}=2 E 0 / E c ( 0 ) = 2 R − 3 ≪ ( E 0 / E c ( 0 ) ) 2 − 1 R^{-3}\ll(E_0/E_{\mathrm{c}}^{(0)})^2-1 R − 3 ≪ ( E 0 / E c ( 0 ) ) 2 − 1 Θ m ≃ 54.7 ∘ \Theta_{\mathrm m}\simeq54.7^\circ Θ m ≃ 54. 7 ∘
At Θ = Θ ∗ \Theta=\Theta_* Θ = Θ ∗
The condition δ / λ 0 ≪ 1 \delta/\lambda_0\ll1 δ / λ 0 ≪ 1 R − 3 ≪ ( E 0 / E c ( 0 ) ) 2 − 1 R^{-3}\ll(E_0/E_c^{(0)})^2-1 R − 3 ≪ ( E 0 / E c ( 0 ) ) 2 − 1 Θ m = cos − 1 ( 1 / 3 ) ≃ 54.7 ∘ \Theta_m=\cos^{-1}(1/\sqrt{3})\simeq54.7^\circ Θ m = cos − 1 ( 1/ 3 ) ≃ 54. 7 ∘ Θ ∗ \Theta_* Θ ∗ δ ≪ λ 0 ≪ 1 \delta\ll\lambda_0\ll1 δ ≪ λ 0 ≪ 1 λ 0 = O ( δ ) \lambda_0=O(\delta) λ 0 = O ( δ )
Long-time behavior of the regular approximation
We first return to the regular expansion with fixed azimuthal angles ϕ 1 \phi_1 ϕ 1 ϕ 2 \phi_2 ϕ 2 P i ( 1 ) P_i^{(1)} P i ( 1 ) (10) , before slow reorientation is absorbed into the leading term as in Eq. (34) . We integrate Eqs. (14) and (15) , form the sum P i ( 0 ) + δ P i ( 1 ) P_i^{(0)}+\delta P_i^{(1)} P i ( 0 ) + δ P i ( 1 ) ϵ 21 = − 0.1097 \epsilon_{21}=-0.1097 ϵ 21 = − 0.1097 σ 21 = − 0.5 \sigma_{21}=-0.5 σ 21 = − 0.5 λ 0 = 1 \lambda_0=1 λ 0 = 1 E 0 / E c ( 0 ) = 2 E_0/E_c^{(0)}=\sqrt{2} E 0 / E c ( 0 ) = 2 R ^ = e x \hat{\mathbf{R}}=\mathbf{e}_x R ^ = e x ϕ 1 = ϕ 2 = π / 6 \phi_1=\phi_2=\pi/6 ϕ 1 = ϕ 2 = π /6 P i ( 1 ) ( 0 ) = 0 P_i^{(1)}(0)=0 P i ( 1 ) ( 0 ) = 0
For these fixed reference phases, Eq. (31) gives a phase change that grows at the slow-time rate
G i , ϕ P ⊥ ( 0 ) = 3 4 sin ( 2 ϕ i ) = 3 3 8 . (48) \frac{G_{i,\phi}}{P_\perp^{(0)}}=\frac{3}{4}\sin(2\phi_i)=\frac{3\sqrt{3}}{8}. \tag*{(48)} P ⊥ ( 0 ) G i , ϕ = 4 3 sin ( 2 ϕ i ) = 8 3 3 . (48) The actual dipole correction grows at rate δ G i , ϕ \delta G_{i,\phi} δ G i , ϕ P i ( 0 ) + δ P i ( 1 ) P_i^{(0)}+\delta P_i^{(1)} P i ( 0 ) + δ P i ( 1 ) Figure 3 (a)]. The unexpanded transverse magnitude approaches a constant value [Figure 3 (a), inset]. In panel (b), the phase change ϕ 1 − ϕ 1 ( 0 ) \phi_1-\phi_1(0) ϕ 1 − ϕ 1 ( 0 ) π / 3 \pi/3 π /3 ϕ 1 ( 0 ) = π / 6 \phi_1(0)=\pi/6 ϕ 1 ( 0 ) = π /6 π / 2 \pi/2 π /2 (37) reduces here to d ϕ / d τ = ( 3 / 4 ) sin ( 2 ϕ ) \mathrm{d}\phi/\mathrm{d}\tau=(3/4)\sin(2\phi) d ϕ / d τ = ( 3/4 ) sin ( 2 ϕ ) Figure 3 (b)].
Figure 3. Loss of uniform validity of the regular expansion at long times for R ^ = e x \widehat{\boldsymbol{R}}=\boldsymbol{e}_x R = e x λ 0 = 1 \lambda_0=1 λ 0 = 1 ϕ 1 ( 0 ) = ϕ 2 ( 0 ) = π / 6 \phi_1(0)=\phi_2(0)=\pi/6 ϕ 1 ( 0 ) = ϕ 2 ( 0 ) = π /6 P ⊥ ( 0 ) + δ P ⊥ ( 1 ) P_{\perp}^{(0)}+\delta P_{\perp}^{(1)} P ⊥ ( 0 ) + δ P ⊥ ( 1 ) P ⊥ ( 0 ) P_{\perp}^{(0)} P ⊥ ( 0 ) δ P 1 , ϕ ( 1 ) / P ⊥ ( 0 ) \delta P_{1,\phi}^{(1)}/P_{\perp}^{(0)} δ P 1 , ϕ ( 1 ) / P ⊥ ( 0 ) π / 3 \pi/3 π /3 π / 2 \pi/2 π /2
Next choose ϕ 1 = ϕ 2 = π / 2 \phi_1=\phi_2=\pi/2 ϕ 1 = ϕ 2 = π /2 G i , ϕ = 0 G_{i,\phi}=0 G i , ϕ = 0 P i , ρ ( 1 ) = 0.1190 P_{i,\rho}^{(1)}=0.1190 P i , ρ ( 1 ) = 0.1190 P i , z ( 1 ) = 0.2593 P_{i,z}^{(1)}=0.2593 P i , z ( 1 ) = 0.2593 Figure 4 ). These steady coefficients follow from Eqs. (29) and (30) . Thus its steady transverse magnitude is P ⊥ ( 0 ) + 0.1190 δ + O ( δ 2 ) P_\perp^{(0)}+0.1190\delta+O(\delta^2) P ⊥ ( 0 ) + 0.1190 δ + O ( δ 2 ) ϵ 21 − M a + 0.2593 δ + O ( δ 2 ) \epsilon_{21}-\mathrm{Ma}+0.2593\delta+O(\delta^2) ϵ 21 − Ma + 0.2593 δ + O ( δ 2 )
Figure 4. Bounded first-order coefficients when the phases satisfy the solvability condition: ϕ 1 = ϕ 2 = π / 2 \phi_1=\phi_2=\pi/2 ϕ 1 = ϕ 2 = π /2 R ^ = e x \widehat{\boldsymbol{R}}=\boldsymbol{e}_x R = e x Figure 3 . Dotted lines are the steady values from Eqs. (29) and (30) ; the phase component vanishes. The dipole correction is δ \delta δ
Failure of the regular expansion near onset
We now examine the steady first-order magnitude correction as the isolated critical field is approached, λ 0 → 0 + \lambda_0\to0^+ λ 0 → 0 + (28) , so that G i , ϕ = 0 G_{i,\phi}=0 G i , ϕ = 0 δ P i , ρ ( 1 ) \delta P_{i,\rho}^{(1)} δ P i , ρ ( 1 ) P ⊥ ( 0 ) P_\perp^{(0)} P ⊥ ( 0 )
We use the first and third component equations in Eq. (27) . Setting their time derivatives to zero gives
0 = − P ⊥ ( 0 ) M a P i , z ( 1 ) + G i , ρ , 0=-\frac{P_\perp^{(0)}}{\mathrm{Ma}}P_{i,z}^{(1)}+G_{i,\rho}, 0 = − Ma P ⊥ ( 0 ) P i , z ( 1 ) + G i , ρ , 0 = 2 P ⊥ ( 0 ) M a P i , ρ ( 1 ) − P i , z ( 1 ) + G i , z . (49) 0=\frac{2P_\perp^{(0)}}{\mathrm{Ma}}P_{i,\rho}^{(1)}-P_{i,z}^{(1)}+G_{i,z}. \tag*{(49)} 0 = Ma 2 P ⊥ ( 0 ) P i , ρ ( 1 ) − P i , z ( 1 ) + G i , z . (49) These are the two linear equations previously solved to obtain Eqs. (29) and (30) . Their coefficient matrix is the ( ρ , z ) (\rho,z) ( ρ , z ) (26) , with determinant
det [ J i ] ( ρ , z ) = 2 ( P ⊥ ( 0 ) ) 2 M a 2 = 2 λ 0 . (50) \det[\mathcal{J}_i]_{(\rho,z)}=\frac{2(P_\perp^{(0)})^2}{\mathrm{Ma}^2}=2\lambda_0. \tag*{(50)} det [ J i ] ( ρ , z ) = Ma 2 2 ( P ⊥ ( 0 ) ) 2 = 2 λ 0 . (50) The homogeneous perturbations governed by this block have growth rates satisfying s 2 + s + 2 λ 0 = 0 s^2+s+2\lambda_0=0 s 2 + s + 2 λ 0 = 0 λ 0 → 0 + \lambda_0\to0^+ λ 0 → 0 + − 1 -1 − 1 − 2 λ 0 -2\lambda_0 − 2 λ 0
To determine the size of that correction, we must also evaluate the forcing. The first equation in Eq. (49) gives
P i , z ( 1 ) = M a P ⊥ ( 0 ) G i , ρ , (51) P_{i,z}^{(1)}=\frac{\mathrm{Ma}}{P_\perp^{(0)}}G_{i,\rho}, \tag*{(51)} P i , z ( 1 ) = P ⊥ ( 0 ) Ma G i , ρ , (51) which is Eq. (29) . Substituting this expression into the second equation and using ( P ⊥ ( 0 ) ) 2 = M a 2 λ 0 (P_\perp^{(0)})^2=\mathrm{Ma}^2\lambda_0 ( P ⊥ ( 0 ) ) 2 = Ma 2 λ 0
2 λ 0 P i , ρ ( 1 ) = G i , ρ − P ⊥ ( 0 ) M a G i , z . (52) 2\lambda_0P_{i,\rho}^{(1)}=G_{i,\rho}-\frac{P_\perp^{(0)}}{\mathrm{Ma}}G_{i,z}. \tag*{(52)} 2 λ 0 P i , ρ ( 1 ) = G i , ρ − Ma P ⊥ ( 0 ) G i , z . (52) This is Eq. (30) rearranged to display the small factor multiplying the unknown magnitude correction.
The right-hand side also becomes small near onset. To evaluate it, substitute the isolated rotating dipoles from Eq. (21) into the forcing in Eq. (25) . Using Eq. (20) to write
M a = M a c 1 + λ 0 , P ⊥ ( 0 ) = M a λ 0 , (53) \mathrm{Ma}=\frac{\mathrm{Ma}_c}{1+\lambda_0},\qquad P_\perp^{(0)}=\mathrm{Ma}\sqrt{\lambda_0}, \tag*{(53)} Ma = 1 + λ 0 Ma c , P ⊥ ( 0 ) = Ma λ 0 , (53) and expanding at fixed pair orientation gives
G i , ρ = P ⊥ ( 0 ) M a G i , z = P ⊥ ( 0 ) [ − 2 σ 21 e z ⋅ Π ⋅ e z − 1 2 e ϕ i ⋅ Π ⋅ e ϕ j ] + O ( ( P ⊥ ( 0 ) ) 2 ) . (54) G_{i,\rho}=\frac{P_{\perp}^{(0)}}{\mathrm{Ma}}G_{i,z}
=P_{\perp}^{(0)}\left[-2\sigma_{21}\boldsymbol{e}_z\cdot\boldsymbol{\Pi}\cdot\boldsymbol{e}_z-\frac{1}{2}\boldsymbol{e}_{\phi_i}\cdot\boldsymbol{\Pi}\cdot\boldsymbol{e}_{\phi_j}\right]+\mathcal{O}\left((P_{\perp}^{(0)})^2\right). \tag*{(54)} G i , ρ = Ma P ⊥ ( 0 ) G i , z = P ⊥ ( 0 ) [ − 2 σ 21 e z ⋅ Π ⋅ e z − 2 1 e ϕ i ⋅ Π ⋅ e ϕ j ] + O ( ( P ⊥ ( 0 ) ) 2 ) . (54) Dividing by 2 λ 0 = 2 ( P ⊥ ( 0 ) ) 2 / M a 2 2\lambda_0=2(P_{\perp}^{(0)})^2/\mathrm{Ma}^2 2 λ 0 = 2 ( P ⊥ ( 0 ) ) 2 / Ma 2
P i , ρ ( 1 ) = M a 2 2 P ⊥ ( 0 ) [ − 2 σ 21 e z ⋅ Π ⋅ e z − 1 2 e ϕ i ⋅ Π ⋅ e ϕ j ] + O ( 1 ) . (55) P_{i,\rho}^{(1)}=\frac{\mathrm{Ma}^2}{2P_{\perp}^{(0)}}\left[-2\sigma_{21}\boldsymbol{e}_z\cdot\boldsymbol{\Pi}\cdot\boldsymbol{e}_z-\frac{1}{2}\boldsymbol{e}_{\phi_i}\cdot\boldsymbol{\Pi}\cdot\boldsymbol{e}_{\phi_j}\right]+\mathcal{O}(1). \tag*{(55)} P i , ρ ( 1 ) = 2 P ⊥ ( 0 ) Ma 2 [ − 2 σ 21 e z ⋅ Π ⋅ e z − 2 1 e ϕ i ⋅ Π ⋅ e ϕ j ] + O ( 1 ) . (55) When the coefficient in brackets approaches a nonzero value, P i , ρ ( 1 ) = O ( λ 0 − 1 / 2 ) P_{i,\rho}^{(1)}=\mathcal{O}(\lambda_0^{-1/2}) P i , ρ ( 1 ) = O ( λ 0 − 1/2 ) (25) gives G i , ρ = O ( P ⊥ ( 0 ) ) G_{i,\rho}=\mathcal{O}(P_{\perp}^{(0)}) G i , ρ = O ( P ⊥ ( 0 ) ) (29) shows that the longitudinal correction P i , z ( 1 ) P_{i,z}^{(1)} P i , z ( 1 )
For these generic configurations, the transverse correction relative to the isolated transverse dipole has size
δ ∣ P i , ρ ( 1 ) ∣ P ⊥ ( 0 ) = O ( δ λ 0 ) . (56) \frac{\delta\lvert P_{i,\rho}^{(1)}\rvert}{P_{\perp}^{(0)}}=\mathcal{O}\left(\frac{\delta}{\lambda_0}\right). \tag*{(56)} P ⊥ ( 0 ) δ ∣ P i , ρ ( 1 ) ∣ = O ( λ 0 δ ) . (56) The regular expansion therefore requires λ 0 ≫ δ \lambda_0\gg\delta λ 0 ≫ δ λ 0 = O ( δ ) \lambda_0=\mathcal{O}(\delta) λ 0 = O ( δ ) P ⊥ ( 0 ) = M a λ 0 P^{(0)}_{\perp}=\mathrm{Ma}\sqrt{\lambda_0} P ⊥ ( 0 ) = Ma λ 0
P ⊥ ( 0 ) = O ( δ 1 / 2 ) . (57) P^{(0)}_{\perp}=\mathcal{O}(\delta^{1/2}). \tag*{(57)} P ⊥ ( 0 ) = O ( δ 1/2 ) . (57) If the coefficient in brackets in Eq. (55) vanishes at onset, the leading singular term is absent and further terms must be examined for that configuration.
In Figure 5 (a), the blue curve is the first-order coefficient from the steady correction equations, and the dashed curve is its λ 0 − 1 / 2 \lambda_0^{-1/2} λ 0 − 1/2 λ 0 \lambda_0 λ 0 Figure 4 , now varying λ 0 \lambda_0 λ 0 e x \boldsymbol{e}_x e x e y \boldsymbol{e}_y e y Π = diag ( − 2 , 1 , 1 ) \boldsymbol{\Pi}=\operatorname{diag}(-2,1,1) Π = diag ( − 2 , 1 , 1 ) e ϕ i = e ϕ j = − e x \boldsymbol{e}_{\phi_i}=\boldsymbol{e}_{\phi_j}=-\boldsymbol{e}_x e ϕ i = e ϕ j = − e x (17) has no e y \boldsymbol{e}_y e y P i , ρ ( 1 ) P^{(1)}_{i,\rho} P i , ρ ( 1 ) σ 21 = − 1 / 2 \sigma_{21}=-1/2 σ 21 = − 1/2 (55) equals 2, confirming that the leading singular term does not cancel. Figure 5 (a) shows the resulting growth of P i , ρ ( 1 ) P^{(1)}_{i,\rho} P i , ρ ( 1 ) λ 0 − 1 / 2 \lambda_0^{-1/2} λ 0 − 1/2
Figure 5. Breakdown of the regular expansion for the transverse dipole magnitude near onset for R ^ = e x \hat{\boldsymbol{R}}=\boldsymbol{e}_x R ^ = e x ϕ 1 = ϕ 2 = π / 2 \phi_1=\phi_2=\pi/2 ϕ 1 = ϕ 2 = π /2 ∣ P 1 , ρ ( 1 ) − P n r , ρ ( 1 ) ∣ \left|P^{(1)}_{1,\rho}-P^{(1)}_{nr,\rho}\right| P 1 , ρ ( 1 ) − P n r , ρ ( 1 ) M a c / λ 0 \mathrm{Ma}_c/\sqrt{\lambda_0} Ma c / λ 0 δ ∣ P 1 , ρ ( 1 ) ∣ / P ⊥ ( 0 ) \delta\left|P^{(1)}_{1,\rho}\right|/P^{(0)}_{\perp} δ P 1 , ρ ( 1 ) / P ⊥ ( 0 ) λ 0 \lambda_0 λ 0 Figure 4 fixed.
The consequence for the regular expansion is shown more directly in Figure 5 (b), which plots the actual correction relative to the leading transverse dipole, δ ∣ P i , ρ ( 1 ) ∣ / P ⊥ ( 0 ) \delta\left|P^{(1)}_{i,\rho}\right|/P^{(0)}_{\perp} δ P i , ρ ( 1 ) / P ⊥ ( 0 ) λ 0 \lambda_0 λ 0 λ 0 ≫ δ \lambda_0\gg\delta λ 0 ≫ δ λ 0 = O ( δ ) \lambda_0=\mathcal{O}(\delta) λ 0 = O ( δ )
On the isolated rotating branch, the transverse dipole magnitude and angular speed vanish as λ 0 → 0 + \lambda_0\to0^+ λ 0 → 0 + σ 21 \sigma_{21} σ 21 P = σ 21 e z \boldsymbol{P}=\sigma_{21}\boldsymbol{e}_z P = σ 21 e z
Near-onset expansion about the critical non-rotating state
Detuning, time scales, and expansion
§ 1.4 identifies the balance λ 0 = O ( δ ) \lambda_0=\mathcal{O}(\delta) λ 0 = O ( δ )
Using the isolated critical Mason number M a c \mathrm{Ma}_c Ma c (3) , introduce a small parameter η \eta η ρ \rho ρ
η = δ 1 / 2 , ρ = λ 0 η 2 , M a − 1 = 1 + ρ η 2 M a c . (58) \eta=\delta^{1/2}, \qquad\rho=\frac{\lambda_0}{\eta^2}, \qquad\mathrm{Ma}^{-1}=\frac{1+\rho\eta^2}{\mathrm{Ma}_c}. \tag*{(58)} η = δ 1/2 , ρ = η 2 λ 0 , Ma − 1 = Ma c 1 + ρ η 2 . (58) Here ρ \rho ρ ρ = [ ( E 0 / E c ) 2 − 1 ] / η 2 \rho=[(E_0/E_c)^2-1]/\eta^2 ρ = [( E 0 / E c ) 2 − 1 ] / η 2 η → 0 \eta\to0 η → 0 ρ \rho ρ ρ \rho ρ [61 ] .
Introduce the time scales
T 0 = t , T 1 = η t , T 2 = η 2 t , T_0=t,\qquad T_1=\eta t,\qquad T_2=\eta^2t, T 0 = t , T 1 = η t , T 2 = η 2 t , d d t = D 0 + η D 1 + η 2 D 2 , D n = ∂ ∂ T n . (59) \frac{\mathrm{d}}{\mathrm{d}t}=D_0+\eta D_1+\eta^2D_2,\qquad D_n=\frac{\partial}{\partial T_n}. \tag*{(59)} d t d = D 0 + η D 1 + η 2 D 2 , D n = ∂ T n ∂ . (59) We initially allow dependence on all three scales. The solvability conditions will determine which slow scales are needed. The expansion for the dipole is
P i = σ 21 e z + η P i [ 1 ] + η 2 P i [ 2 ] + η 3 P i [ 3 ] + O ( η 4 ) . (60) P_i=\sigma_{21}e_z+\eta P_i^{[1]}+\eta^2P_i^{[2]}+\eta^3P_i^{[3]}+O(\eta^4). \tag*{(60)} P i = σ 21 e z + η P i [ 1 ] + η 2 P i [ 2 ] + η 3 P i [ 3 ] + O ( η 4 ) . (60) where each coefficient may depend on ( T 0 , T 1 , T 2 ) (T_0,T_1,T_2) ( T 0 , T 1 , T 2 ) η \eta η δ \delta δ (10) . In particular, P i [ 2 ] P_i^{[2]} P i [ 2 ] η 2 = δ \eta^2=\delta η 2 = δ P i [ 1 ] P_i^{[1]} P i [ 1 ] η \eta η (57) .
We substitute Eqs. (58) –(60) into each of Eqs. (7). The local-field equation immediately gives
E i = e z − η 2 σ 21 Π ⋅ e z − η 3 Π ⋅ P j [ 1 ] + O ( η 4 ) . (61) E_i=e_z-\eta^2\sigma_{21}\boldsymbol{\Pi}\cdot e_z-\eta^3\boldsymbol{\Pi}\cdot P_j^{[1]}+O(\eta^4). \tag*{(61)} E i = e z − η 2 σ 21 Π ⋅ e z − η 3 Π ⋅ P j [ 1 ] + O ( η 4 ) . (61) We also expand the angular velocity as Ω i = η ω i [ 1 ] + η 2 ω i [ 2 ] + η 3 ω i [ 3 ] + ⋯ \Omega_i=\eta\omega_i^{[1]}+\eta^2\omega_i^{[2]}+\eta^3\omega_i^{[3]}+\cdots Ω i = η ω i [ 1 ] + η 2 ω i [ 2 ] + η 3 ω i [ 3 ] + ⋯ η \eta η
Order eta: the leading transverse dipole
The torque and dipole equations give
ω i [ 1 ] = Ma c − 1 P i [ 1 ] × e z , D 0 P i [ 1 ] = − ( P i [ 1 ] ⋅ e z ) e z . (62) \omega_i^{[1]}=\operatorname{Ma}_c^{-1}P_i^{[1]}\times e_z,\qquad D_0P_i^{[1]}=-(P_i^{[1]}\cdot e_z)e_z. \tag*{(62)} ω i [ 1 ] = Ma c − 1 P i [ 1 ] × e z , D 0 P i [ 1 ] = − ( P i [ 1 ] ⋅ e z ) e z . (62) Hence
P i [ 1 ] = a i ( T 1 , T 2 ) + c i ( T 1 , T 2 ) e − T 0 e z , a i ⋅ e z = 0. (63) P_i^{[1]}=a_i(T_1,T_2)+c_i(T_1,T_2)e^{-T_0}e_z,\qquad a_i\cdot e_z=0. \tag*{(63)} P i [ 1 ] = a i ( T 1 , T 2 ) + c i ( T 1 , T 2 ) e − T 0 e z , a i ⋅ e z = 0. (63) The longitudinal perturbation decays on the charge-relaxation time. After this initial transient, P i [ 1 ] = a i ( T 1 , T 2 ) P_i^{[1]}=a_i(T_1,T_2) P i [ 1 ] = a i ( T 1 , T 2 ) η a i \eta a_i η a i a 1 a_1 a 1 a 2 a_2 a 2
Order eta squared: the first solvability condition
At this order the torque balance gives
ω i [ 2 ] = Ma c − 1 [ P i [ 2 ] + σ 21 2 Π ⋅ e z ] × e z . (64) \omega_i^{[2]}=\operatorname{Ma}_c^{-1}\left[P_i^{[2]}+\sigma_{21}^2\boldsymbol{\Pi}\cdot e_z\right]\times e_z. \tag*{(64)} ω i [ 2 ] = Ma c − 1 [ P i [ 2 ] + σ 21 2 Π ⋅ e z ] × e z . (64) Using this expression in the dipole equation gives
D 0 P i [ 2 ] + D 1 a i = [ ∣ a i ∣ 2 Ma c − P i [ 2 ] ] ⋅ e z − σ 21 2 e z ⋅ Π ⋅ e z . (65) D_0P_i^{[2]}+D_1a_i=
\left[\frac{|a_i|^2}{\operatorname{Ma}_c}-P_i^{[2]}\right]\cdot e_z
-\sigma_{21}^2e_z\cdot\boldsymbol{\Pi}\cdot e_z. \tag*{(65)} D 0 P i [ 2 ] + D 1 a i = [ Ma c ∣ a i ∣ 2 − P i [ 2 ] ] ⋅ e z − σ 21 2 e z ⋅ Π ⋅ e z . (65) The right-hand side has no transverse component. Thus
D 0 P i , ⊥ [ 2 ] = − D 1 a i . (66) D_0P_{i,\perp}^{[2]}=-D_1a_i. \tag*{(66)} D 0 P i , ⊥ [ 2 ] = − D 1 a i . (66) Its solution contains − T 0 D 1 a i -T_0D_1a_i − T 0 D 1 a i
D 1 a i = 0 , a i = a i ( T 2 ) . (67) D_1a_i=0,\qquad a_i=a_i(T_2). \tag*{(67)} D 1 a i = 0 , a i = a i ( T 2 ) . (67) The leading transverse amplitude a i a_i a i T 1 = η t T_1=\eta t T 1 = η t T 2 T_2 T 2
The longitudinal part of Eq. (65) relaxes to a value determined by ∣ a i ∣ 2 |a_i|^2 ∣ a i ∣ 2 (66) determines the transverse part of P i [ 2 ] P_i^{[2]} P i [ 2 ] T 0 T_0 T 0 η a i \eta a_i η a i η 2 \eta^2 η 2 η \eta η a i a_i a i a i a_i a i
P i [ 2 ] = − σ 21 2 Π ⋅ e z + ∣ a i ∣ 2 Ma c e z , ω i [ 2 ] = 0. (68) P_i^{[2]}=-\sigma_{21}^2\boldsymbol{\Pi}\cdot e_z+\frac{|a_i|^2}{\operatorname{Ma}_c}e_z,\qquad\omega_i^{[2]}=0. \tag*{(68)} P i [ 2 ] = − σ 21 2 Π ⋅ e z + Ma c ∣ a i ∣ 2 e z , ω i [ 2 ] = 0. (68) The first term is the electric interaction correction to the non-rotating dipole found in Eq. (17) , with δ = η 2 \delta=\eta^2 δ = η 2 a i = 0 a_i=0 a i = 0
Order eta cubed: evolution on the slow time
The torque coefficient at this order is
ω i [ 3 ] = Ma c − 1 [ P i [ 3 ] × e z + ρ a i × e z − σ 21 a i × ( Π ⋅ e z ) − σ 21 e z × ( Π ⋅ a j ) − 1 2 Π ⋅ ( a j × e z ) ] . (69) \omega_i^{[3]}=\operatorname{Ma}_c^{-1}\left[P_i^{[3]}\times e_z+\rho a_i\times e_z-\sigma_{21}a_i\times(\boldsymbol{\Pi}\cdot e_z)-\sigma_{21}e_z\times(\boldsymbol{\Pi}\cdot a_j)-\frac{1}{2}\boldsymbol{\Pi}\cdot(a_j\times e_z)\right]. \tag*{(69)} ω i [ 3 ] = Ma c − 1 [ P i [ 3 ] × e z + ρ a i × e z − σ 21 a i × ( Π ⋅ e z ) − σ 21 e z × ( Π ⋅ a j ) − 2 1 Π ⋅ ( a j × e z ) ] . (69) Substitute this expression and Eq. (68) into the order-η 3 \eta^3 η 3
D 0 P i , ⊥ [ 3 ] = F i − D 2 a i , (70) D_0P_{i,\perp}^{[3]}=\mathcal{F}_i-D_2a_i, \tag*{(70)} D 0 P i , ⊥ [ 3 ] = F i − D 2 a i , (70) where the known forcing is
F i = [ ρ − 2 σ 21 ( e z ⋅ Π ⋅ e z ) ] a i − ∣ a i ∣ 2 Ma c 2 a i + 1 2 e z × [ Π ⋅ ( e z × a i ) ] . (71) \mathcal{F}_i=\left[\rho-2\sigma_{21}(e_z\cdot\boldsymbol{\Pi}\cdot e_z)\right]a_i-\frac{|a_i|^2}{\operatorname{Ma}_c^2}a_i+\frac{1}{2}e_z\times\left[\boldsymbol{\Pi}\cdot(e_z\times a_i)\right]. \tag*{(71)} F i = [ ρ − 2 σ 21 ( e z ⋅ Π ⋅ e z ) ] a i − Ma c 2 ∣ a i ∣ 2 a i + 2 1 e z × [ Π ⋅ ( e z × a i ) ] . (71) The right-hand side of Eq. (70) is independent of T 0 T_0 T 0 F i − D 2 a i \mathcal{F}_i-D_2\mathbf{a}_i F i − D 2 a i T 0 T_0 T 0 P i , ⊥ [ 3 ] P_{i,\perp}^{[3]} P i , ⊥ [ 3 ] D 2 a i = F i D_2\mathbf{a}_i=\mathcal{F}_i D 2 a i = F i
d a i d T 2 = [ ρ − 2 σ 21 ( e z ⋅ Π ⋅ e z ) ] a i − ∣ a i ∣ 2 M a c 2 a i + 1 2 e z × [ Π ⋅ ( e z × a j ) ] . (72) \frac{d\mathbf{a}_i}{dT_2}=\left[\rho-2\sigma_{21}(\mathbf{e}_z\cdot\mathbf{\Pi}\cdot\mathbf{e}_z)\right]\mathbf{a}_i-\frac{|\mathbf{a}_i|^2}{\mathrm{Ma}_c^2}\mathbf{a}_i+\frac{1}{2}\mathbf{e}_z\times\left[\mathbf{\Pi}\cdot(\mathbf{e}_z\times\mathbf{a}_j)\right]. \tag*{(72)} d T 2 d a i = [ ρ − 2 σ 21 ( e z ⋅ Π ⋅ e z ) ] a i − Ma c 2 ∣ a i ∣ 2 a i + 2 1 e z × [ Π ⋅ ( e z × a j ) ] . (72) Equation (72) is the central result of the near-onset analysis. It determines both the magnitude and direction of each leading transverse dipole, η a i \eta\mathbf{a}_i η a i T 2 = η 2 t = δ t T_2=\eta^2t=\delta t T 2 = η 2 t = δ t ρ a i \rho\mathbf{a}_i ρ a i − 2 σ 21 ( e z ⋅ Π ⋅ e z ) a i -2\sigma_{21}(\mathbf{e}_z\cdot\mathbf{\Pi}\cdot\mathbf{e}_z)\mathbf{a}_i − 2 σ 21 ( e z ⋅ Π ⋅ e z ) a i
The negative cubic term comes from the longitudinal dipole correction in Eq. (68) , generated by convection of surface charge by each particle’s own rotation. This local nonlinear feedback limits the growth of the transverse dipole and is already present for an isolated rotating sphere. Like the detuning and electric interaction terms, it acts parallel to a i \mathbf{a}_i a i
The final term is the hydrodynamic interaction with the other particle. It depends on a j \mathbf{a}_j a j Π \mathbf{\Pi} Π a i \mathbf{a}_i a i
Each a i \mathbf{a}_i a i (72) therefore consists of four coupled scalar ordinary differential equations, one for each of X 1 X_1 X 1 Y 1 Y_1 Y 1 X 2 X_2 X 2 Y 2 Y_2 Y 2 a i = X i e x + Y i e y \mathbf{a}_i=X_i\mathbf{e}_x+Y_i\mathbf{e}_y a i = X i e x + Y i e y
The local-field derivative in Eq. (7c) contributes at order η 3 \eta^3 η 3 D 0 P j [ 1 ] D_0P_j^{[1]} D 0 P j [ 1 ] (67) gives ϵ 21 E ˙ i = − ϵ 21 η 5 Π ⋅ d a j / d T 2 + ⋯ \epsilon_{21}\dot{\mathbf{E}}_i=-\epsilon_{21}\eta^5\mathbf{\Pi}\cdot d\mathbf{a}_j/dT_2+\cdots ϵ 21 E ˙ i = − ϵ 21 η 5 Π ⋅ d a j / d T 2 + ⋯ (72) .
Connection with the phase dynamics away from onset
For a nonzero amplitude, write a i = ∣ a i ∣ e ρ i \mathbf{a}_i=|\mathbf{a}_i|\mathbf{e}_{\rho_i} a i = ∣ a i ∣ e ρ i (72) along e ϕ i \mathbf{e}_{\phi_i} e ϕ i
∣ a i ∣ d ϕ i d T 2 = ∣ a j ∣ 2 e ρ i ⋅ Π ⋅ e ϕ j . (73) |\mathbf{a}_i|\frac{d\phi_i}{dT_2}=\frac{|\mathbf{a}_j|}{2}\mathbf{e}_{\rho_i}\cdot\mathbf{\Pi}\cdot\mathbf{e}_{\phi_j}. \tag*{(73)} ∣ a i ∣ d T 2 d ϕ i = 2 ∣ a j ∣ e ρ i ⋅ Π ⋅ e ϕ j . (73) For equal leading magnitudes, the constant-phase condition agrees with the P ⊥ ( 0 ) → 0 + P_{\perp}^{(0)}\to0^+ P ⊥ ( 0 ) → 0 + (28) . The near-onset equations thus recover the limiting phase condition while also determining the magnitude, which can no longer be prescribed by the isolated-sphere solution.
Without interactions, Eq. (72) reduces to d a i / d T 2 = ρ a i − ∣ a i ∣ 2 a i / M a c 2 d\mathbf{a}_i/dT_2=\rho\mathbf{a}_i-|\mathbf{a}_i|^2\mathbf{a}_i/\mathrm{Ma}_c^2 d a i / d T 2 = ρ a i − ∣ a i ∣ 2 a i / Ma c 2 ρ > 0 \rho>0 ρ > 0
∣ P i , ⊥ ∣ = M a c λ 0 + O ( λ 0 3 / 2 ) , (74) |P_{i,\perp}|=\mathrm{Ma}_c\sqrt{\lambda_0}+O(\lambda_0^{3/2}), \tag*{(74)} ∣ P i , ⊥ ∣ = Ma c λ 0 + O ( λ 0 3/2 ) , (74) in agreement with P ⊥ ( 0 ) = M a λ 0 P_{\perp}^{(0)}=\mathrm{Ma}\sqrt{\lambda_0} P ⊥ ( 0 ) = Ma λ 0
M a = M a c 1 + λ 0 , M a λ 0 = M a c λ 0 + O ( λ 0 3 / 2 ) . (75) \mathrm{Ma}=\frac{\mathrm{Ma}_c}{1+\lambda_0},\qquad\mathrm{Ma}\sqrt{\lambda_0}=\mathrm{Ma}_c\sqrt{\lambda_0}+O(\lambda_0^{3/2}). \tag*{(75)} Ma = 1 + λ 0 Ma c , Ma λ 0 = Ma c λ 0 + O ( λ 0 3/2 ) . (75) Four rotation patterns near onset
Linear growth rates and the four modes
Choose the line of centres in the x x x z z z R ^ = sin Θ e x + cos Θ e z \hat{\mathbf{R}}=\sin\Theta\mathbf{e}_x+\cos\Theta\mathbf{e}_z R ^ = sin Θ e x + cos Θ e z Θ ≠ 0 \Theta\ne0 Θ = 0 e x \mathbf{e}_x e x e y \mathbf{e}_y e y a i = X i e x + Y i e y \mathbf{a}_i=X_i\mathbf{e}_x+Y_i\mathbf{e}_y a i = X i e x + Y i e y (72) gives
d X i d T 2 = g X i − 1 2 X j − X i 2 + Y i 2 M a c 2 X i , (76) \frac{dX_i}{dT_2}=gX_i-\frac{1}{2}X_j-\frac{X_i^2+Y_i^2}{\mathrm{Ma}_c^2}X_i, \tag*{(76)} d T 2 d X i = g X i − 2 1 X j − Ma c 2 X i 2 + Y i 2 X i , (76) d Y i d T 2 = g Y i + h Y j − X i 2 + Y i 2 M a c 2 Y i , (77) \frac{dY_i}{dT_2}=gY_i+hY_j-\frac{X_i^2+Y_i^2}{\mathrm{Ma}_c^2}Y_i, \tag*{(77)} d T 2 d Y i = g Y i + h Y j − Ma c 2 X i 2 + Y i 2 Y i , (77) where j = 3 − i j=3-i j = 3 − i
g = ρ − 2 σ 21 ( 1 − 3 cos 2 Θ ) , h = 1 2 ( 2 − 3 cos 2 Θ ) . (78) g=\rho-2\sigma_{21}(1-3\cos^2\Theta),\qquad h=\frac{1}{2}(2-3\cos^2\Theta). \tag*{(78)} g = ρ − 2 σ 21 ( 1 − 3 cos 2 Θ ) , h = 2 1 ( 2 − 3 cos 2 Θ ) . (78) Here g g g − 1 / 2 -1/2 − 1/2 h h h (72) . To find the onset of rotation, linearize about X 1 = Y 1 = X 2 = Y 2 = 0 X_1=Y_1=X_2=Y_2=0 X 1 = Y 1 = X 2 = Y 2 = 0
d d T 2 ( X 1 X 2 ) = ( g − 1 / 2 − 1 / 2 g ) ( X 1 X 2 ) , (79) \frac{d}{dT_2}\begin{pmatrix}X_1\\X_2\end{pmatrix}=
\begin{pmatrix}g&-1/2\\-1/2&g\end{pmatrix}
\begin{pmatrix}X_1\\X_2\end{pmatrix}, \tag*{(79)} d T 2 d ( X 1 X 2 ) = ( g − 1/2 − 1/2 g ) ( X 1 X 2 ) , (79) d d T 2 ( Y 1 Y 2 ) = ( g h h g ) ( Y 1 Y 2 ) . \frac{d}{dT_2}\begin{pmatrix}Y_1\\Y_2\end{pmatrix}=
\begin{pmatrix}g&h\\h&g\end{pmatrix}
\begin{pmatrix}Y_1\\Y_2\end{pmatrix}. d T 2 d ( Y 1 Y 2 ) = ( g h h g ) ( Y 1 Y 2 ) . Exchanging the two identical particles leaves these linearized equations unchanged. Perturbations that are equal on the two particles and those that are opposite therefore evolve independently, as expected from particle-exchange symmetry [63 ] . Adding and subtracting the equations implements this decomposition. Define
X ± = X 1 ± X 2 2 , Y ± = Y 1 ± Y 2 2 . (80) X_{\pm} = \frac{X_1 \pm X_2}{\sqrt{2}}, \qquad Y_{\pm} = \frac{Y_1 \pm Y_2}{\sqrt{2}}. \tag*{(80)} X ± = 2 X 1 ± X 2 , Y ± = 2 Y 1 ± Y 2 . (80) Each coordinate obeys d Z k / d T 2 = Λ k Z k dZ_k/dT_2 = \Lambda_k Z_k d Z k / d T 2 = Λ k Z k Z k = X + , X − , Y + , Y − Z_k = X_+,X_-,Y_+,Y_- Z k = X + , X − , Y + , Y − Λ k = ρ + m k \Lambda_k = \rho+ m_k Λ k = ρ + m k
m X ± = − 2 σ 21 ( 1 − 3 cos 2 Θ ) ∓ 1 2 , (81) m_{X_{\pm}} = -2\sigma_{21}(1-3\cos^2\Theta) \mp\frac{1}{2}, \tag*{(81)} m X ± = − 2 σ 21 ( 1 − 3 cos 2 Θ ) ∓ 2 1 , (81) m Y ± = − 2 σ 21 ( 1 − 3 cos 2 Θ ) ± 1 2 ( 2 − 3 cos 2 Θ ) . (82) m_{Y_{\pm}} = -2\sigma_{21}(1-3\cos^2\Theta) \pm\frac{1}{2}(2-3\cos^2\Theta). \tag*{(82)} m Y ± = − 2 σ 21 ( 1 − 3 cos 2 Θ ) ± 2 1 ( 2 − 3 cos 2 Θ ) . (82) The upper signs give X + X_+ X + Y + Y_+ Y + X − X_- X − Y − Y_- Y − Z k ( T 2 ) = Z k ( 0 ) e Λ k T 2 Z_k(T_2) = Z_k(0)e^{\Lambda_kT_2} Z k ( T 2 ) = Z k ( 0 ) e Λ k T 2 Λ k < 0 \Lambda_k < 0 Λ k < 0 Λ k > 0 \Lambda_k > 0 Λ k > 0
The + sign describes equal transverse dipoles on the two particles; the − sign describes opposite transverse dipoles. After the initial transient, Eq. (62) gives ω i [ 1 ] = M a c − 1 a i × e z \omega_i^{[1]} = \mathrm{Ma}_c^{-1}\boldsymbol{a}_i \times\boldsymbol{e}_z ω i [ 1 ] = Ma c − 1 a i × e z (68) gives ω i [ 2 ] = 0 \omega_i^{[2]} = 0 ω i [ 2 ] = 0 a i = X i e x + Y i e y \boldsymbol{a}_i = X_i\boldsymbol{e}_x + Y_i\boldsymbol{e}_y a i = X i e x + Y i e y
Ω i = η M a c a i × e z + O ( η 3 ) = η M a c ( Y i e x − X i e y ) + O ( η 3 ) . (83) \begin{aligned}
\boldsymbol{\Omega}_i = \frac{\eta}{\mathrm{Ma}_c}\boldsymbol{a}_i \times\boldsymbol{e}_z + O(\eta^3) \\
&= \frac{\eta}{\mathrm{Ma}_c}(Y_i\boldsymbol{e}_x-X_i\boldsymbol{e}_y)+O(\eta^3). \tag*{(83)}
\end{aligned} Ω i = Ma c η a i × e z + O ( η 3 ) = Ma c η ( Y i e x − X i e y ) + O ( η 3 ) . (83) Thus X + X_+ X + Y + Y_+ Y + X − X_- X − Y − Y_- Y − X X X e y \boldsymbol{e}_y e y Y Y Y e x \boldsymbol{e}_x e x Figure 6 illustrates these four patterns.
Figure 6. The four leading rotation patterns near onset. The page is the x x x z z z ± e x \pm e_x ± e x X X X ± e x \pm e_x ± e x Y Y Y ± e y \pm e_y ± e y a k a_k a k (87) .
Setting Λ k = 0 \Lambda_k = 0 Λ k = 0 ρ c ( k ) = − m k \rho_c^{(k)} = -m_k ρ c ( k ) = − m k Figure 7 (a) shows how the four interaction shifts vary with pair angle. The largest shift is m X − m_{X_-} m X − Θ m \Theta_m Θ m m Y + m_{Y_+} m Y + λ 0 = ρ η 2 \lambda_0 = \rho\eta^2 λ 0 = ρ η 2
Figure 7. Growth-rate shifts and steady amplitudes near onset for ϵ 21 = − 0.1097 \epsilon_{21}=-0.1097 ϵ 21 = − 0.1097 σ 21 = − 0.5 \sigma_{21}=-0.5 σ 21 = − 0.5 m k m_k m k Θ m = cos − 1 ( 1 / 3 ) \Theta_m=\cos^{-1}(1/\sqrt{3}) Θ m = cos − 1 ( 1/ 3 ) m X − m_{X_-} m X − m Y + m_{Y_+} m Y + ∣ a k ∣ / M a c |a_k|/M_{\mathrm{ac}} ∣ a k ∣/ M ac Θ = 30 ∘ \Theta=30^\circ Θ = 3 0 ∘
( E c , k E c ( 0 ) ) 2 = 1 − η 2 m k + O ( η 4 ) . (84) \left(\frac{E_{c,k}}{E_c^{(0)}}\right)^2 = 1-\eta^2m_k+O(\eta^4). \tag*{(84)} ( E c ( 0 ) E c , k ) 2 = 1 − η 2 m k + O ( η 4 ) . (84) Equation (84) is the leading near-onset asymptotic prediction for the interaction-induced change in each modal critical field. A positive m k m_k m k m k m_k m k
Since T 2 = η 2 t T_2 = \eta^2t T 2 = η 2 t
e Λ k T 2 = e η 2 Λ k t = e η 2 Λ k t p h y s / τ M W . (85) e^{\Lambda_kT_2}=e^{\eta^2\Lambda_kt}=e^{\eta^2\Lambda_kt_{\mathrm{phys}}/\tau_{\mathrm{MW}}}. \tag*{(85)} e Λ k T 2 = e η 2 Λ k t = e η 2 Λ k t phys / τ MW . (85) Its leading growth rate is therefore η 2 Λ k \eta^2\Lambda_k η 2 Λ k η 2 Λ k / τ M W \eta^2\Lambda_k/\tau_{\mathrm{MW}} η 2 Λ k / τ MW
When Θ = 0 \Theta= 0 Θ = 0 m X ± = m Y ± = 4 σ 21 ∓ 1 / 2 m_{X_{\pm}} = m_{Y_{\pm}} = 4\sigma_{21}\mp1/2 m X ± = m Y ± = 4 σ 21 ∓ 1/2 X X X Y Y Y
Steady amplitudes and their stability
For an X X X X 1 = a k X_1 = a_k X 1 = a k X 2 = q a k X_2 = qa_k X 2 = q a k Y 1 = Y 2 = 0 Y_1 = Y_2 = 0 Y 1 = Y 2 = 0 q = + 1 q = +1 q = + 1 q = − 1 q = -1 q = − 1 q q q X q X_q X q Y q Y_q Y q Y Y Y Y 1 = a k Y_1 = a_k Y 1 = a k Y 2 = q a k Y_2 = qa_k Y 2 = q a k X 1 = X 2 = 0 X_1 = X_2 = 0 X 1 = X 2 = 0 (76) and (77) gives the scalar equation
d a k d T 2 = Λ k a k − a k 3 M a c 2 . (86) \frac{d a_k}{dT_2}=\Lambda_k a_k-\frac{a_k^3}{\mathrm{Ma}_c^2}. \tag*{(86)} d T 2 d a k = Λ k a k − Ma c 2 a k 3 . (86) Besides a k = 0 a_k = 0 a k = 0
a k = ± M a c Λ k , Λ k > 0. (87) a_k=\pm\mathrm{Ma}_c\sqrt{\Lambda_k}, \qquad\Lambda_k>0. \tag*{(87)} a k = ± Ma c Λ k , Λ k > 0. (87) Their magnitude grows continuously from zero as the field exceeds the critical field for that mode. The two signs reverse both particles’ rotations. Here a k a_k a k (80) has magnitude 2 ∣ a k ∣ \sqrt{2}|a_k| 2 ∣ a k ∣ η ∣ a k ∣ \eta|a_k| η ∣ a k ∣ η Λ k \eta\sqrt{\Lambda_k} η Λ k
Equation (86) describes perturbations that keep the chosen rotation pattern. To establish stability we must also perturb the other dipole components. For example, consider an X − X_- X − ( X 1 , X 2 ) = ( a k , − a k ) (X_1,X_2)=(a_k,-a_k) ( X 1 , X 2 ) = ( a k , − a k ) Δ X ( T 2 ) \Delta X(T_2) Δ X ( T 2 ) Δ Y ( T 2 ) \Delta Y(T_2) Δ Y ( T 2 )
perturbation growth rate on T 2 ( Δ X 1 , Δ X 2 ) = ( Δ X , − Δ X ) − 2 Λ X − ( Δ X 1 , Δ X 2 ) = ( Δ X , Δ X ) Λ X + − 3 Λ X − ( Δ Y 1 , Δ Y 2 ) = ( Δ Y , Δ Y ) Λ Y + − Λ X − ( Δ Y 1 , Δ Y 2 ) = ( Δ Y , − Δ Y ) Λ Y − − Λ X − (88) \begin{array}{c|c}
\text{perturbation} & \text{growth rate on }T_2\\
(\Delta X_1,\Delta X_2)=(\Delta X,-\Delta X) & -2\Lambda_{X_-}\\
(\Delta X_1,\Delta X_2)=(\Delta X,\Delta X) & \Lambda_{X_+}-3\Lambda_{X_-}\\
(\Delta Y_1,\Delta Y_2)=(\Delta Y,\Delta Y) & \Lambda_{Y_+}-\Lambda_{X_-}\\
(\Delta Y_1,\Delta Y_2)=(\Delta Y,-\Delta Y) & \Lambda_{Y_-}-\Lambda_{X_-} \tag*{(88)}
\end{array} perturbation ( Δ X 1 , Δ X 2 ) = ( Δ X , − Δ X ) ( Δ X 1 , Δ X 2 ) = ( Δ X , Δ X ) ( Δ Y 1 , Δ Y 2 ) = ( Δ Y , Δ Y ) ( Δ Y 1 , Δ Y 2 ) = ( Δ Y , − Δ Y ) growth rate on T 2 − 2 Λ X − Λ X + − 3 Λ X − Λ Y + − Λ X − Λ Y − − Λ X − (88) These rates follow by substituting the steady state plus the indicated small perturbations into Eqs. (76) and (77) , and retaining terms linear in Δ X \Delta X Δ X Δ Y \Delta Y Δ Y x x x y y y d ( Δ Y ) / d T 2 = ( Λ Y − − Λ X − ) Δ Y d(\Delta Y)/dT_2=(\Lambda_{Y_-}-\Lambda_{X_-})\Delta Y d ( Δ Y ) / d T 2 = ( Λ Y − − Λ X − ) Δ Y
A steady state is stable if all four growth rates are negative. The full set of growth rates for all four patterns is given in Table S1 in Sec. SIII of the SM. It shows that the mode with the unique largest m k m_k m k ρ + m k > 0 \rho+m_k>0 ρ + m k > 0 Θ = 30 ∘ \Theta=30^\circ Θ = 3 0 ∘ X − X_- X − Figure 7 (b) shows that its nonzero branch is stable above this field (solid curve). The other single-mode branches are unstable (dashed curves). Each branch emerges continuously from zero at its own critical detuning and then grows with the square root of the distance above that threshold. The horizontal-axis intercepts therefore give the critical detuning of each branch; solid and dashed lines denote stable and unstable nonzero branches, respectively.
Table S1. Growth rates about each nonzero single-mode rotating branch, with Λ k = ρ + m k > 0 \Lambda_k=\rho+m_k>0 Λ k = ρ + m k > 0 X − X_- X − (88) of the main text.
To determine which pattern first appears, substitute Eqs. (81) and (82) and cancel their common electric contribution:
m X − − m X + = 1 , m Y + − m X − = 1 2 ( 1 − 3 cos 2 Θ ) . (89) m_{X_-}-m_{X_+}=1,\qquad m_{Y_+}-m_{X_-}=\frac{1}{2}(1-3\cos^2\Theta). \tag*{(89)} m X − − m X + = 1 , m Y + − m X − = 2 1 ( 1 − 3 cos 2 Θ ) . (89) Also m Y − ≤ m X − m_{Y_-}\le m_{X_-} m Y − ≤ m X − Θ = 0 \Theta=0 Θ = 0 X − X_- X − 0 < Θ < Θ m 0<\Theta<\Theta_m 0 < Θ < Θ m Y + Y_+ Y + Θ m < Θ ≤ π / 2 \Theta_m<\Theta\le\pi/2 Θ m < Θ ≤ π /2
Θ m = cos − 1 ( 1 / 3 ) ≃ 54.7 ∘ . (90) \Theta_m=\cos^{-1}(1/\sqrt{3})\simeq54.7^\circ. \tag*{(90)} Θ m = cos − 1 ( 1/ 3 ) ≃ 54. 7 ∘ . (90) This angle also agrees with the small-amplitude limit of the away-from-onset phase-selection boundary in Eq. (40) : cos 2 Θ ∗ = 1 / [ 3 ( 1 + λ 0 ) ] \cos^2\Theta_* = 1/[3(1+\lambda_0)] cos 2 Θ ∗ = 1/ [ 3 ( 1 + λ 0 )] 1 / 3 1/3 1/3 λ 0 → 0 + \lambda_0\to0^+ λ 0 → 0 + δ ≪ λ 0 ≪ 1 \delta\ll\lambda_0\ll1 δ ≪ λ 0 ≪ 1
At Θ = 0 \Theta= 0 Θ = 0 Θ = Θ m \Theta= \Theta_m Θ = Θ m X − X_- X − Y + Y_+ Y +
Separating electric and hydrodynamic interactions
Governing equations and non-rotating state
To separate the effects of the two interaction mechanisms, introduce s E , s H ∈ { 0 , 1 } s_{\mathrm{E}},s_{\mathrm{H}} \in\{0,1\} s E , s H ∈ { 0 , 1 }
E i = e z − s E δ Π ⋅ P j , (91) \mathbf{E}_i = \mathbf{e}_z - s_{\mathrm{E}}\delta\mathbf{\Pi}\cdot\mathbf{P}_j, \tag*{(91)} E i = e z − s E δ Π ⋅ P j , (91) Ω i + s H δ 2 Π ⋅ Ω j = M a − 1 P i × E i , (92) \mathbf{\Omega}_i + s_{\mathrm{H}}\frac{\delta}{2}\mathbf{\Pi}\cdot\mathbf{\Omega}_j = \mathrm{Ma}^{-1}\mathbf{P}_i \times\mathbf{E}_i, \tag*{(92)} Ω i + s H 2 δ Π ⋅ Ω j = Ma − 1 P i × E i , (92) P ˙ i + s E ϵ 21 δ Π ⋅ P ˙ j = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) . (93) \dot{\mathbf{P}}_i + s_{\mathrm{E}}\epsilon_{21}\delta\mathbf{\Pi}\cdot\dot{\mathbf{P}}_j = \mathbf{\Omega}_i \times\left(\mathbf{P}_i-\epsilon_{21}\mathbf{E}_i\right)-\left(\mathbf{P}_i-\sigma_{21}\mathbf{E}_i\right). \tag*{(93)} P ˙ i + s E ϵ 21 δ Π ⋅ P ˙ j = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) . (93) The choices ( s E , s H ) = ( 1 , 1 ) , ( 1 , 0 ) , ( 0 , 1 ) , ( 0 , 0 ) (s_{\mathrm{E}},s_{\mathrm{H}})=(1,1),(1,0),(0,1),(0,0) ( s E , s H ) = ( 1 , 1 ) , ( 1 , 0 ) , ( 0 , 1 ) , ( 0 , 0 )
The common non-rotating field and dipole are
E n r = ( I + s E σ 21 δ Π ) − 1 ⋅ e z , P n r = σ 21 E n r . (94) \mathbf{E}_{\mathrm{nr}}=\left(\mathbf{I}+s_{\mathrm{E}}\sigma_{21}\delta\mathbf{\Pi}\right)^{-1}\cdot\mathbf{e}_z,\qquad\mathbf{P}_{\mathrm{nr}}=\sigma_{21}\mathbf{E}_{\mathrm{nr}}. \tag*{(94)} E nr = ( I + s E σ 21 δ Π ) − 1 ⋅ e z , P nr = σ 21 E nr . (94) Because Ω i = 0 \mathbf{\Omega}_i=0 Ω i = 0 s E s_{\mathrm{E}} s E s E = 0 s_{\mathrm{E}}=0 s E = 0 E n r = e z \mathbf{E}_{\mathrm{nr}}=\mathbf{e}_z E nr = e z P n r = σ 21 e z \mathbf{P}_{\mathrm{nr}}=\sigma_{21}\mathbf{e}_z P nr = σ 21 e z
Stationary onset in the unexpanded model
We first obtain the stationary critical fields of the unexpanded model, then compare their first-order expansion with the near-onset asymptotic predictions. In the linear stability analysis of the non-rotating state, write P i = P n r + p i \mathbf{P}_i=\mathbf{P}_{\mathrm{nr}}+\mathbf{p}_i P i = P nr + p i Ω i = ω i \mathbf{\Omega}_i=\boldsymbol{\omega}_i Ω i = ω i (91) then gives E i = E n r − s E δ Π ⋅ p j \mathbf{E}_i=\mathbf{E}_{\mathrm{nr}}-s_{\mathrm{E}}\delta\mathbf{\Pi}\cdot\mathbf{p}_j E i = E nr − s E δ Π ⋅ p j p 2 = q p 1 \mathbf{p}_2=q\mathbf{p}_1 p 2 = q p 1 ω 2 = q ω 1 \boldsymbol{\omega}_2=q\boldsymbol{\omega}_1 ω 2 = q ω 1 q = + 1 q=+1 q = + 1 − 1 -1 − 1 ω q = ω 1 \boldsymbol{\omega}_q=\boldsymbol{\omega}_1 ω q = ω 1
At a stationary critical field the leading growth rate reaches zero. The time derivatives in the linearized Eq. (89c) then vanish, giving p i + s E σ 21 δ Π ⋅ p j = M a c E n r × ω i \mathbf{p}_i+s_{\mathrm{E}}\sigma_{21}\delta\mathbf{\Pi}\cdot\mathbf{p}_j=\mathrm{Ma}_c\mathbf{E}_{\mathrm{nr}}\times\boldsymbol{\omega}_i p i + s E σ 21 δ Π ⋅ p j = Ma c E nr × ω i M a − 1 = ( E 0 / E c ( 0 ) ) 2 / M a c \mathrm{Ma}^{-1}=(E_0/E_c^{(0)})^2/\mathrm{Ma}_c Ma − 1 = ( E 0 / E c ( 0 ) ) 2 / Ma c (20) , gives the neutral angular-velocity equation
( I + q s H δ 2 Π ) ⋅ ω q = ( E 0 E c ( 0 ) ) 2 [ ∣ E n r ∣ 2 I − E n r E n r ⋅ ] ω q . (95) \left(\mathbf{I}+\frac{q s_{\mathrm{H}}\delta}{2}\mathbf{\Pi}\right)\cdot\boldsymbol{\omega}_q=\left(\frac{E_0}{E_c^{(0)}}\right)^2\left[\lvert\mathbf{E}_{\mathrm{nr}}\rvert^2\mathbf{I}-\mathbf{E}_{\mathrm{nr}}\mathbf{E}_{\mathrm{nr}}\cdot\right]\boldsymbol{\omega}_q. \tag*{(95)} ( I + 2 q s H δ Π ) ⋅ ω q = ( E c ( 0 ) E 0 ) 2 [ ∣ E nr ∣ 2 I − E nr E nr ⋅ ] ω q . (95) For each sign of q q q δ \delta δ
Because E n r \mathbf{E}_{\mathrm{nr}} E nr R ^ \hat{\mathbf{R}} R ^ x x x z z z X q X_q X q X + X_+ X + q = + 1 q=+1 q = + 1 X − X_- X − q = − 1 q=-1 q = − 1 Y q Y_q Y q X q X_q X q Y q Y_q Y q (95) and gives
( E c , X q E c ( 0 ) ) 2 = 1 + q s H δ / 2 ∣ E n r ∣ 2 , (96) \left(\frac{E_{\mathrm{c},X_q}}{E_c^{(0)}}\right)^2=\frac{1+q s_{\mathrm{H}}\delta/2}{\lvert\mathbf{E}_{\mathrm{nr}}\rvert^2}, \tag*{(96)} ( E c ( 0 ) E c , X q ) 2 = ∣ E nr ∣ 2 1 + q s H δ /2 , (96) ( E c , Y q E c ( 0 ) ) 2 = ( 1 + q s H δ / 2 ) ( 1 − q s H δ ) ( 1 + q s H δ / 2 ) ∣ E n r ∣ 2 − ( 3 q s H δ / 2 ) ( R ^ ⋅ E n r ) 2 . (97) \left(\frac{E_{\mathrm{c},Y_q}}{E_c^{(0)}}\right)^2=\frac{(1+q s_{\mathrm{H}}\delta/2)(1-q s_{\mathrm{H}}\delta)}{(1+q s_{\mathrm{H}}\delta/2)\lvert\mathbf{E}_{\mathrm{nr}}\rvert^2-(3q s_{\mathrm{H}}\delta/2)(\hat{\mathbf{R}}\cdot\mathbf{E}_{\mathrm{nr}})^2}. \tag*{(97)} ( E c ( 0 ) E c , Y q ) 2 = ( 1 + q s H δ /2 ) ∣ E nr ∣ 2 − ( 3 q s H δ /2 ) ( R ^ ⋅ E nr ) 2 ( 1 + q s H δ /2 ) ( 1 − q s H δ ) . (97) These expressions are exact within the unexpanded dipole–rotlet model. The first onset occurs at the smallest of the four positive modal fields. The local-field derivative term in Eq. (89c) affects nonzero growth rates but vanishes in the stationary neutral condition. Normalization by E c ( 0 ) E_c^{(0)} E c ( 0 ) ϵ 21 \epsilon_{21} ϵ 21 (96) and (97) , so these normalized stationary thresholds depend on the conductivity contrast through E n r \mathbf{E}_{\mathrm{nr}} E nr s E s_{\mathrm{E}} s E s H s_{\mathrm{H}} s H
m X ± = − 2 s E σ 21 ( 1 − 3 cos 2 Θ ) ± s H 2 , (98) m_{X_\pm}=-2s_{\mathrm{E}}\sigma_{21}(1-3\cos^2\Theta)\pm\frac{s_{\mathrm{H}}}{2}, \tag*{(98)} m X ± = − 2 s E σ 21 ( 1 − 3 cos 2 Θ ) ± 2 s H , (98) m Y ± = − 2 s E σ 21 ( 1 − 3 cos 2 Θ ) ± s H 2 ( 2 − 3 cos 2 Θ ) . (99) m_{Y_\pm}=-2s_{\mathrm{E}}\sigma_{21}(1-3\cos^2\Theta)\pm\frac{s_{\mathrm{H}}}{2}(2-3\cos^2\Theta). \tag*{(99)} m Y ± = − 2 s E σ 21 ( 1 − 3 cos 2 Θ ) ± 2 s H ( 2 − 3 cos 2 Θ ) . (99) For comparison with the unexpanded model, recall δ = η 2 \delta=\eta^2 δ = η 2 (96) and (97) to first order in δ \delta δ m k m_k m k (84) gives each modal squared critical field. The lowest of these fields corresponds to the largest m k m_k m k
E c E c ( 0 ) = 1 − δ max k m k + O ( δ 2 ) . (100) \frac{E_c}{E_c^{(0)}}=\sqrt{1-\delta\max_k m_k+O(\delta^2)}. \tag*{(100)} E c ( 0 ) E c = 1 − δ k max m k + O ( δ 2 ) . (100) Equation (100) therefore contains the same first-order asymptotic information as Eq. (84) , with the switches included and the first mode selected. Figure 8 (a) separates the two mechanisms: electric interaction alone raises onset near a field-aligned pair and lowers it near a perpendicular pair for this material, while hydrodynamic coupling splits the competing rotation modes. Their combined effect depends on pair angle. Panel (b) resolves the four modal fields and identifies which one sets the first-onset envelope. The solid and dashed curves remain close because the leading asymptotic shifts are common to both calculations; their differences arise from higher powers of δ \delta δ
Figure 8. Stationary onset fields of the dipole–rotlet model at R = 2.5 R = 2.5 R = 2.5 ε 21 = − 0.1097 \varepsilon_{21} = -0.1097 ε 21 = − 0.1097 σ 21 = − 0.5 \sigma_{21} = -0.5 σ 21 = − 0.5 (100) . (b) Four modal thresholds with both interactions. The dotted black curve gives the lowest of the four thresholds at each angle. The vertical dotted line marks the unexpanded-model X − − Y + X_--Y_+ X − − Y + Θ c r o s s = 50.8 ∘ \Theta_{\mathrm{cross}} = 50.8^\circ Θ cross = 50. 8 ∘
Interaction mechanisms and finite-separation crossover
We now use the unexpanded thresholds (92)–(93) to distinguish the two interaction mechanisms and determine how finite separation changes the crossover between the first unstable modes. With electric interaction only, s H = 0 s_H=0 s H = 0 E c / E c ( 0 ) = ∣ E n r ∣ − 1 E_c/E_c^{(0)}=|\mathbf{E}_{\mathrm{nr}}|^{-1} E c / E c ( 0 ) = ∣ E nr ∣ − 1 E n r = e z \mathbf{E}_{\mathrm{nr}}=\mathbf{e}_z E nr = e z X q X_q X q 1 + q δ / 2 1+q\delta/2 1 + q δ /2 Y q Y_q Y q (97) with ∣ E n r ∣ = 1 |\mathbf{E}_{\mathrm{nr}}|=1 ∣ E nr ∣ = 1 R ^ ⋅ E n r = cos Θ \hat{\mathbf{R}}\cdot\mathbf{E}_{\mathrm{nr}}=\cos\Theta R ^ ⋅ E nr = cos Θ X − X_- X − Y + Y_+ Y +
When the hydrodynamic interaction is present, the first instability is controlled by X − X_- X − Y + Y_+ Y + X + X_+ X + Y − Y_- Y − X − X_- X − X − X_- X − Y + Y_+ Y + Θ c r o s s ( R ) \Theta_{\mathrm{cross}}(R) Θ cross ( R ) X − X_- X − Y + Y_+ Y +
( R ^ ⋅ E n r ) 2 ∣ E n r ∣ 2 = 1 + δ / 2 3 ( 1 − δ / 2 ) . (101) \frac{\left(\hat{\mathbf{R}}\cdot\mathbf{E}_{\mathrm{nr}}\right)^2}{|\mathbf{E}_{\mathrm{nr}}|^2}
=\frac{1+\delta/2}{3(1-\delta/2)}. \tag*{(101)} ∣ E nr ∣ 2 ( R ^ ⋅ E nr ) 2 = 3 ( 1 − δ /2 ) 1 + δ /2 . (101) In the leading near-onset asymptotic limit, this reduces to cos 2 Θ = 1 / 3 \cos^2\Theta=1/3 cos 2 Θ = 1/3 δ \delta δ
cos 2 Θ c r o s s H ( R ) = 1 + 1 / ( 2 R 3 ) 3 [ 1 − 1 / ( 2 R 3 ) ] . (102) \cos^2\Theta_{\mathrm{cross}}^{\mathrm{H}}(R)
=\frac{1+1/(2R^3)}{3[1-1/(2R^3)]}. \tag*{(102)} cos 2 Θ cross H ( R ) = 3 [ 1 − 1/ ( 2 R 3 )] 1 + 1/ ( 2 R 3 ) . (102) The left-hand side of Eq. (101) is the squared cosine of the angle between the line of centres and the local non-rotating electric field. With hydrodynamic interaction alone, this field is the imposed field, so the left-hand side is cos 2 Θ \cos^2\Theta cos 2 Θ Θ \Theta Θ σ 21 = − 0.5 \sigma_{21}=-0.5 σ 21 = − 0.5 R = 2.5 R=2.5 R = 2.5 Θ c r o s s = 50.8 ∘ \Theta_{\mathrm{cross}}=50.8^\circ Θ cross = 50. 8 ∘ Θ c r o s s H = 53.4 ∘ \Theta_{\mathrm{cross}}^{\mathrm{H}}=53.4^\circ Θ cross H = 53. 4 ∘ Θ m = 54.7 ∘ \Theta_m=54.7^\circ Θ m = 54. 7 ∘ R = 6 R=6 R = 6 54.45 ∘ 54.45^\circ 54.4 5 ∘ 54.64 ∘ 54.64^\circ 54.6 4 ∘ Figure 8 (a,Figure 8 b) results from taking the lower of the smooth X − X_- X − Y + Y_+ Y +
Steady rotation away from onset
The near-onset amplitude equations apply when λ 0 = O ( δ ) \lambda_0=\mathcal{O}(\delta) λ 0 = O ( δ ) λ 0 > 0 \lambda_0>0 λ 0 > 0 δ / λ 0 ≪ 1 \delta/\lambda_0\ll1 δ / λ 0 ≪ 1 E 0 / E c ( 0 ) = 2 E_0/E_c^{(0)}=2 E 0 / E c ( 0 ) = 2 (40) is Θ ∗ ≃ 73.2 ∘ \Theta_* \simeq73.2^\circ Θ ∗ ≃ 73. 2 ∘ (42) is stable above it. We therefore evaluate the counter-rotating state ( ϕ 1 , ϕ 2 ) = ( 0 , π ) (\phi_1,\phi_2)=(0,\pi) ( ϕ 1 , ϕ 2 ) = ( 0 , π ) Θ = 0 , π / 4 \Theta=0,\pi/4 Θ = 0 , π /4 ϕ 1 = ϕ 2 = π / 2 \phi_1=\phi_2=\pi/2 ϕ 1 = ϕ 2 = π /2 Θ = π / 2 \Theta=\pi/2 Θ = π /2 X − X_- X − Y + Y_+ Y +
For the theory and simulation comparisons, define the mean angular speed
Ω ‾ = ∣ Ω 1 ∣ + ∣ Ω 2 ∣ 2 . (103) \overline{\Omega}=\frac{|\Omega_1|+|\Omega_2|}{2}. \tag*{(103)} Ω = 2 ∣ Ω 1 ∣ + ∣ Ω 2 ∣ . (103) We impose Eq. (28) , substitute the steady dipole corrections (29) and (30) into Eq. (13) , and obtain the first-order regular asymptotic predictions
Ω ‾ X = λ 0 − δ ( P ⊥ ( 0 ) ) 2 + M a 2 P ⊥ ( 0 ) M a [ ( ϵ 21 − M a ) ( 1 − 3 cos 2 Θ ) − 1 4 ] + O ( δ 2 ) . (104) \overline{\Omega}_X
=\sqrt{\lambda_0}
-\delta\frac{(P_{\perp}^{(0)})^2+\mathrm{Ma}^2}{P_{\perp}^{(0)}\mathrm{Ma}}
\left[(\epsilon_{21}-\mathrm{Ma})(1-3\cos^2\Theta)-\frac{1}{4}\right]
+O(\delta^2). \tag*{(104)} Ω X = λ 0 − δ P ⊥ ( 0 ) Ma ( P ⊥ ( 0 ) ) 2 + Ma 2 [ ( ϵ 21 − Ma ) ( 1 − 3 cos 2 Θ ) − 4 1 ] + O ( δ 2 ) . (104) Ω ‾ Y + ∣ Θ = π / 2 = λ 0 − δ ( P ⊥ ( 0 ) ) 2 + M a 2 P ⊥ ( 0 ) M a ( ε 21 − M a − 1 2 ) + O ( δ 2 ) . (105) \left.\overline{\Omega}_{Y_+}\right|_{\Theta=\pi/2}
= \sqrt{\lambda_0}
- \delta\frac{\left(P_{\perp}^{(0)}\right)^2+\mathrm{Ma}^2}{P_{\perp}^{(0)}\mathrm{Ma}}
\left(\varepsilon_{21}-\mathrm{Ma}-\frac{1}{2}\right)
+ O\left(\delta^2\right). \tag*{(105)} Ω Y + Θ = π /2 = λ 0 − δ P ⊥ ( 0 ) Ma ( P ⊥ ( 0 ) ) 2 + Ma 2 ( ε 21 − Ma − 2 1 ) + O ( δ 2 ) . (105) In each bracket, the term containing ε 21 − M a \varepsilon_{21}-\mathrm{Ma} ε 21 − Ma − 1 / 4 -1/4 − 1/4 − 1 / 2 -1/2 − 1/2
Boundary element formulation
The boundary element method [64 , 65 ] is well suited to simulating the electrohydrodynamics of particles. It solves Laplace’s equation for the electric potential [58 , 66 ] and the Stokes equations for the fluid flow [67 , 68 ] through boundary integral equations, requiring only the particle surfaces to be discretized. Here, we extend the validated single-particle solver [48 ] to two interacting rigid spheres. The simulations evolve the surface charge and resolve the coupled electric and hydrodynamic interactions, providing a numerical benchmark for the dipole–rotlet model and asymptotic theory. For the Stokes problem, we use the completed, deflated double-layer formulation of Pozrikidis [69 ] . Figure 9 shows the flat triangular mesh used to discretize each sphere.
Figure 9. Discretized sphere pair at R = 2.5 R = 2.5 R = 2.5 Θ = π / 4 \Theta= \pi/4 Θ = π /4
To extend the validated single-particle solver to two spheres with minimal changes, we retain its electroviscous time scale τ E = μ / ( ε 1 E 0 2 ) \tau_E = \mu/(\varepsilon_1 E_0^2) τ E = μ / ( ε 1 E 0 2 ) t p h y s t_{\mathrm{phys}} t phys Ω p h y s \boldsymbol{\Omega}_{\mathrm{phys}} Ω phys t B = t p h y s / τ E t_B = t_{\mathrm{phys}}/\tau_E t B = t phys / τ E Ω B = τ E Ω p h y s \boldsymbol{\Omega}_B = \tau_E\boldsymbol{\Omega}_{\mathrm{phys}} Ω B = τ E Ω phys M a = 2 τ E / τ M W \mathrm{Ma}=2\tau_E/\tau_{\mathrm{MW}} Ma = 2 τ E / τ MW
t = M a 2 t B , Ω = 2 M a Ω B . (106) t=\frac{\mathrm{Ma}}{2}t_B,\qquad\Omega=\frac{2}{\mathrm{Ma}}\Omega_B. \tag*{(106)} t = 2 Ma t B , Ω = Ma 2 Ω B . (106) Let S = S 1 ∪ S 2 S=S_1\cup S_2 S = S 1 ∪ S 2 s i ( x ) = x − x c , i \boldsymbol{s}_i(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{x}_{c,i} s i ( x ) = x − x c , i i i i a a a E 0 E_0 E 0 ϵ 1 E 0 \epsilon_1E_0 ϵ 1 E 0
Q = ϵ 2 ϵ 1 , κ = σ 2 σ 1 , χ = Q + 2 κ + 2 . (107) Q=\frac{\epsilon_2}{\epsilon_1},\qquad\kappa=\frac{\sigma_2}{\sigma_1},\qquad\chi=\frac{Q+2}{\kappa+2}. \tag*{(107)} Q = ϵ 1 ϵ 2 , κ = σ 1 σ 2 , χ = κ + 2 Q + 2 . (107) Superscripts + + + − - − E t + = E t − \boldsymbol{E}_t^+=\boldsymbol{E}_t^- E t + = E t − q s = E n + − Q E n − q_s=E_n^+-QE_n^- q s = E n + − Q E n −
Δ E n , i = E n , i + − E n , i − . (108) \Delta E_{n,i}=E_{n,i}^+-E_{n,i}^-. \tag*{(108)} Δ E n , i = E n , i + − E n , i − . (108) It is also the density of the single-layer potential. With G ( x 0 , x ) = [ 4 π ∣ x 0 − x ∣ ] − 1 G(\boldsymbol{x}_0,\boldsymbol{x})=[4\pi|\boldsymbol{x}_0-\boldsymbol{x}|]^{-1} G ( x 0 , x ) = [ 4 π ∣ x 0 − x ∣ ] − 1
Φ ( x 0 ) = − e z ⋅ x 0 + ∑ j = 1 2 ∫ S j G ( x 0 , x ) Δ E n , j ( x ) d S . (109) \Phi(\boldsymbol{x}_0)=-\boldsymbol{e}_z\cdot\boldsymbol{x}_0+\sum_{j=1}^{2}\int_{S_j}G(\boldsymbol{x}_0,\boldsymbol{x})\Delta E_{n,j}(\boldsymbol{x})\,\mathrm{d}S. \tag*{(109)} Φ ( x 0 ) = − e z ⋅ x 0 + j = 1 ∑ 2 ∫ S j G ( x 0 , x ) Δ E n , j ( x ) d S . (109) For a target x 0 ∈ S i \boldsymbol{x}_0\in S_i x 0 ∈ S i
∫ S i [ Δ E n , i ( x ) − Δ E n , i ( x 0 ) ] n i ( x 0 ) ⋅ ∇ 0 G d S + ∑ j ≠ i ∫ S j Δ E n , j ( x ) n i ( x 0 ) ⋅ ∇ 0 G d S + Δ E n , i ( x 0 ) [ Q Q − 1 − L i ( x 0 ) ] = e z ⋅ n i ( x 0 ) + q s , i ( x 0 ) Q − 1 . (110) \begin{aligned}
&\int_{S_i}[\Delta E_{n,i}(\boldsymbol{x})-\Delta E_{n,i}(\boldsymbol{x}_0)]\boldsymbol{n}_i(\boldsymbol{x}_0)\cdot\nabla_0G\,\mathrm{d}S\\
&+\sum_{j\ne i}\int_{S_j}\Delta E_{n,j}(\boldsymbol{x})\boldsymbol{n}_i(\boldsymbol{x}_0)\cdot\nabla_0G\,\mathrm{d}S\\
&+\Delta E_{n,i}(\boldsymbol{x}_0)\left[\frac{Q}{Q-1}-L_i(\boldsymbol{x}_0)\right]\\
&=\boldsymbol{e}_z\cdot\boldsymbol{n}_i(\boldsymbol{x}_0)+\frac{q_{s,i}(\boldsymbol{x}_0)}{Q-1}. \tag*{(110)}
\end{aligned} ∫ S i [ Δ E n , i ( x ) − Δ E n , i ( x 0 )] n i ( x 0 ) ⋅ ∇ 0 G d S + j = i ∑ ∫ S j Δ E n , j ( x ) n i ( x 0 ) ⋅ ∇ 0 G d S + Δ E n , i ( x 0 ) [ Q − 1 Q − L i ( x 0 ) ] = e z ⋅ n i ( x 0 ) + Q − 1 q s , i ( x 0 ) . (110) Subtracting the target value of the density regularizes the self-surface integral. The integral over the other sphere has no coincident-point singularity. The geometric correction from the self-surface subtraction of Sellier [70 ] and Das and Saintillan [58 ] takes the following form on sphere i i i
L i ( x 0 ) = n i ( x 0 ) ⋅ ∫ S i { 2 G n i ( x ) + [ n i ( x ) − n i ( x 0 ) ] [ n i ( x ) ⋅ ∇ x G ] } d S . (111) L_i(\boldsymbol{x}_0)=\boldsymbol{n}_i(\boldsymbol{x}_0)\cdot\int_{S_i}\left\{2G\boldsymbol{n}_i(\boldsymbol{x})+[\boldsymbol{n}_i(\boldsymbol{x})-\boldsymbol{n}_i(\boldsymbol{x}_0)][\boldsymbol{n}_i(\boldsymbol{x})\cdot\nabla_xG]\right\}\,\mathrm{d}S. \tag*{(111)} L i ( x 0 ) = n i ( x 0 ) ⋅ ∫ S i { 2 G n i ( x ) + [ n i ( x ) − n i ( x 0 )] [ n i ( x ) ⋅ ∇ x G ] } d S . (111) Here ∇ x G \nabla_xG ∇ x G x \boldsymbol{x} x
After solving Eq. (110) , the normal fields follow from
E n , i + = q s , i − Q Δ E n , i 1 − Q , E n , i − = q s , i − Δ E n , i 1 − Q . (112) E_{n,i}^+=\frac{q_{s,i}-Q\Delta E_{n,i}}{1-Q},\qquad E_{n,i}^-=\frac{q_{s,i}-\Delta E_{n,i}}{1-Q}. \tag*{(112)} E n , i + = 1 − Q q s , i − Q Δ E n , i , E n , i − = 1 − Q q s , i − Δ E n , i . (112) The surface potential is then evaluated using Eq. (109) , and its surface gradient gives
E t , i = − ∇ s Φ . (113) \boldsymbol{E}_{t,i}=-\nabla_s\Phi. \tag*{(113)} E t , i = − ∇ s Φ. (113) This avoids direct evaluation of the strongly singular tangential-field integral. The reconstructed fields determine the dimensionless Maxwell-stress jump and electric torque:
f E , i = [ E i + E i + − 1 2 ∣ E i + ∣ 2 I − Q ( E i − E i − − 1 2 ∣ E i − ∣ 2 I ) ] ⋅ n i , (114) \boldsymbol{f}_{E,i}=\left[\boldsymbol{E}_i^+\boldsymbol{E}_i^+-\frac{1}{2}|\boldsymbol{E}_i^+|^2\boldsymbol{I}-Q\left(\boldsymbol{E}_i^-\boldsymbol{E}_i^--\frac{1}{2}|\boldsymbol{E}_i^-|^2\boldsymbol{I}\right)\right]\cdot\boldsymbol{n}_i, \tag*{(114)} f E , i = [ E i + E i + − 2 1 ∣ E i + ∣ 2 I − Q ( E i − E i − − 2 1 ∣ E i − ∣ 2 I ) ] ⋅ n i , (114) T E , i = ∫ S i s i ( x ) × f E , i d S . (115) \boldsymbol{T}_{E,i}=\int_{S_i}\boldsymbol{s}_i(\boldsymbol{x})\times\boldsymbol{f}_{E,i}\,\mathrm{d}S. \tag*{(115)} T E , i = ∫ S i s i ( x ) × f E , i d S . (115) The steady non-rotating charge distribution gives zero net electric torque on each sphere. Numerical integration on the mesh leaves a small residual T E , i n r \boldsymbol{T}_{E,i}^{\mathrm{nr}} T E , i nr
The exterior fluid obeys the incompressible Stokes equations and is at rest at infinity. To isolate torque-driven rotational coupling, we prescribe the electric torque on each sphere and set both force inputs to zero in the Stokes mobility problem. Let ψ j \boldsymbol{\psi}_j ψ j S j S_j S j [69 , 71 ] gives the exterior velocity at x 0 \boldsymbol{x}_0 x 0
u ( x 0 ) = ∑ j = 1 2 ∫ S j K j ⋅ ψ j ( x ) d S + v ( x 0 ) . (116) \boldsymbol{u}(\boldsymbol{x}_0)=\sum_{j=1}^{2}\int_{S_j}\mathbf{K}_j\cdot\boldsymbol{\psi}_j(\boldsymbol{x})\,\mathrm{d}S+\boldsymbol{v}(\boldsymbol{x}_0). \tag*{(116)} u ( x 0 ) = j = 1 ∑ 2 ∫ S j K j ⋅ ψ j ( x ) d S + v ( x 0 ) . (116) At the integration point x ∈ S j \boldsymbol{x}\in S_j x ∈ S j r = x − x 0 \boldsymbol{r}=\boldsymbol{x}-\boldsymbol{x}_0 r = x − x 0 n j ( x ) \boldsymbol{n}_j(\boldsymbol{x}) n j ( x )
K j ( x , x 0 ) = − 6 r ⋅ n j ( x ) ∣ r ∣ 5 r r . (117) \mathbf{K}_j(\boldsymbol{x},\boldsymbol{x}_0)=-6\frac{\boldsymbol{r}\cdot\boldsymbol{n}_j(\boldsymbol{x})}{|\boldsymbol{r}|^5}\boldsymbol{r}\boldsymbol{r}. \tag*{(117)} K j ( x , x 0 ) = − 6 ∣ r ∣ 5 r ⋅ n j ( x ) r r . (117) The point torques at the sphere centres give the completion velocity
v ( x 0 ) = 1 4 π M a ∑ j = 1 2 ( T E , j − T E , j n r ) × s j ( x 0 ) ∣ s j ( x 0 ) ∣ 3 . (118) \boldsymbol{v}(\boldsymbol{x}_0)=\frac{1}{4\pi\mathrm{Ma}}\sum_{j=1}^{2}\frac{(\boldsymbol{T}_{E,j}-\boldsymbol{T}_{E,j}^{\mathrm{nr}})\times\boldsymbol{s}_j(\boldsymbol{x}_0)}{|\boldsymbol{s}_j(\boldsymbol{x}_0)|^3}. \tag*{(118)} v ( x 0 ) = 4 π Ma 1 j = 1 ∑ 2 ∣ s j ( x 0 ) ∣ 3 ( T E , j − T E , j nr ) × s j ( x 0 ) . (118) The prefactor uses the Maxwell–Wagner time scale and the Mason number in Eq. (2) , which is twice that used in Ref. [48 ] .
At a target x 0 ∈ S i \boldsymbol{x}_0\in S_i x 0 ∈ S i
u ( x 0 ) = U i + Ω i × s i ( x 0 ) . (119) \boldsymbol{u}(\boldsymbol{x}_0)=\boldsymbol{U}_i+\boldsymbol{\Omega}_i\times\boldsymbol{s}_i(\boldsymbol{x}_0). \tag*{(119)} u ( x 0 ) = U i + Ω i × s i ( x 0 ) . (119) Let A i A_i A i S i S_i S i
Pozrikidis [69 ] gives
− ∫ S i K i ⋅ [ ψ i ( x ) − ψ i ( x 0 ) ] d S − ∑ j ≠ i ∫ S j K j ⋅ ψ j ( x ) d S + 4 π A i n i ( x 0 ) ∫ S i ψ i ⋅ n i d S + U i + Ω i × s i ( x 0 ) = v ( x 0 ) . (120) -\int_{S_i} \mathbf{K}_i\cdot[\psi_i(\mathbf{x})-\psi_i(\mathbf{x}_0)]\,\mathrm{d}S
-\sum_{j\ne i}\int_{S_j}\mathbf{K}_j\cdot\psi_j(\mathbf{x})\,\mathrm{d}S
+\frac{4\pi}{A_i}\mathbf{n}_i(\mathbf{x}_0)\int_{S_i}\psi_i\cdot\mathbf{n}_i\,\mathrm{d}S
+\mathbf{U}_i+\boldsymbol{\Omega}_i\times\mathbf{s}_i(\mathbf{x}_0)=\mathbf{v}(\mathbf{x}_0). \tag*{(120)} − ∫ S i K i ⋅ [ ψ i ( x ) − ψ i ( x 0 )] d S − j = i ∑ ∫ S j K j ⋅ ψ j ( x ) d S + A i 4 π n i ( x 0 ) ∫ S i ψ i ⋅ n i d S + U i + Ω i × s i ( x 0 ) = v ( x 0 ) . (120) The local term − 4 π ψ i ( x 0 ) -4\pi\psi_i(\mathbf{x}_0) − 4 π ψ i ( x 0 ) n i ( x 0 ) \mathbf{n}_i(\mathbf{x}_0) n i ( x 0 ) S i S_i S i
The translational and angular velocities in Eq. (113b) follow from the surface moments of the same density [69 ] ,
U i = − 4 π A i ∫ S i ψ i d S , (121) \mathbf{U}_i=-\frac{4\pi}{A_i}\int_{S_i}\psi_i\,\mathrm{d}S, \tag*{(121)} U i = − A i 4 π ∫ S i ψ i d S , (121) I i ⋅ Ω i = − 4 π ∫ S i s i ( x ) × ψ i d S , (122) \mathbf{I}_i\cdot\boldsymbol{\Omega}_i=-4\pi\int_{S_i}\mathbf{s}_i(\mathbf{x})\times\psi_i\,\mathrm{d}S, \tag*{(122)} I i ⋅ Ω i = − 4 π ∫ S i s i ( x ) × ψ i d S , (122) I i = ∫ S i ( ∣ s i ( x ) ∣ 2 I − s i ( x ) s i ( x ) ) d S . (123) \mathbf{I}_i=\int_{S_i}\left(\lvert\mathbf{s}_i(\mathbf{x})\rvert^2\mathbf{I}-\mathbf{s}_i(\mathbf{x})\mathbf{s}_i(\mathbf{x})\right)\,\mathrm{d}S. \tag*{(123)} I i = ∫ S i ( ∣ s i ( x ) ∣ 2 I − s i ( x ) s i ( x ) ) d S . (123) Here I i \mathbf{I}_i I i
D q s , i D t = χ ( κ E n , i − − E n , i + ) , (124) \frac{Dq_{s,i}}{Dt}=\chi(\kappa E_{n,i}^{-}-E_{n,i}^{+}), \tag*{(124)} D t D q s , i = χ ( κ E n , i − − E n , i + ) , (124) d x d t = Ω i × s i ( x ) , x ∈ S i . (125) \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\boldsymbol{\Omega}_i\times\mathbf{s}_i(\mathbf{x}),\qquad\mathbf{x}\in S_i. \tag*{(125)} d t d x = Ω i × s i ( x ) , x ∈ S i . (125) The derivative follows material surface nodes. Rigid rotation preserves area, so their motion incorporates charge convection.
The surfaces use flat P1 triangles, with continuous piecewise-linear densities and 642 nodes and 1280 triangles on each sphere for the main results. Singular element integrals use Duffy-transformed quadrature [72 ] , and the linear systems are solved with GMRES [73 ] . Tangential fields are obtained by local surface differentiation of the potential. Charge and orientation are advanced together by explicit midpoint stepping with Δ t B = 0.1 \Delta t_B=0.1 Δ t B = 0.1
P i = 1 4 π ∫ S i Δ E n , i ( x ) s i ( x ) d S . (126) \mathbf{P}_i=\frac{1}{4\pi}\int_{S_i}\Delta E_{n,i}(\mathbf{x})\mathbf{s}_i(\mathbf{x})\,\mathrm{d}S. \tag*{(126)} P i = 4 π 1 ∫ S i Δ E n , i ( x ) s i ( x ) d S . (126) For the counter-rotating state at R = 8 R=8 R = 8 Θ = 0 \Theta=0 Θ = 0 E 0 / E c ( 0 ) = 2 E_0/E_c^{(0)}=2 E 0 / E c ( 0 ) = 2 (128) , increasing the mesh from 642 to 2562 nodes per sphere changed the steady mean angular speed by 0.035% and the interaction-induced speed correction by approximately 2%. Halving Δ t B \Delta t_B Δ t B 7 × 10 − 11 7\times10^{-11} 7 × 1 0 − 11 R = 2.5 R=2.5 R = 2.5 Θ = 0 \Theta=0 Θ = 0 π / 4 \pi/4 π /4 π / 2 \pi/2 π /2 (127) .
Numerical results
Onset envelope from perturbation growth and decay
The onset calculations use the material ratios Q ≃ 1.4216 Q\simeq1.4216 Q ≃ 1.4216 κ = 10 − 6 \kappa=10^{-6} κ = 1 0 − 6
ϵ 21 ≃ 0.1232 , σ 21 ≃ − 0.5. (127) \epsilon_{21}\simeq0.1232,\qquad\sigma_{21}\simeq-0.5. \tag*{(127)} ϵ 21 ≃ 0.1232 , σ 21 ≃ − 0.5. (127) This is a nearly insulating particle with higher permittivity than the suspending fluid. For each comparison, the BEM simulations, unexpanded dipole–rotlet model and asymptotic predictions use matching material parameters. The angular-speed comparisons below use the second material in Eq. (128) . The small finite conductivity approximates an insulating sphere in both sets of simulations.
We locate onset by following small perturbations of the interacting non-rotating charge distribution. For each geometry, a lower test field is placed 3% below the smaller of the unexpanded-model and near-onset asymptotic first-onset fields, and an upper test field 3% above the larger of the two. All four rotation patterns decay at the lower field, while at least one grows at the upper field. The two fields therefore define an interval containing the onset of rotation. These tests provide an initial comparison with the predicted critical fields. Refining the interval by time integration becomes costly near onset, because the leading growth or decay rate approaches zero and an appreciable change takes progressively longer to observe. A steady-state BEM solver with pseudo-arclength continuation and bifurcation detection could locate onset directly and follow both stable and unstable rotating branches. Each continuation step would require a nonlinear solve coupling surface charge and rigid-body rotation, making an efficient implementation a substantial computational undertaking.
Figure 10 compares these intervals with the near-onset asymptotic and unexpanded-model onset curves at R = 2.5 R=2.5 R = 2.5 E 0 / E c ( 0 ) = 1 E_0/E_c^{(0)}=1 E 0 / E c ( 0 ) = 1 R = 2.5 R = 2.5 R = 2.5 Θ = 0 \Theta= 0 Θ = 0 R = 2.5 R = 2.5 R = 2.5 R = 3 R = 3 R = 3 Θ = π / 2 \Theta= \pi/2 Θ = π /2 R = 4 R = 4 R = 4
Figure 10. Finite-resolution onset range with boundary element simulations for the material ratios in Eq. (127) . Solid curves are the stationary thresholds of the unexpanded model; dashed curves are the leading near-onset asymptotic prediction. Downward and upward open triangles denote BEM fields at which all four tested rotation patterns decay or at least one grows, respectively. BEM simulations are shown at 9° angular intervals.
The near-onset asymptotic curve and the stationary onset curve of the unexpanded dipole–rotlet model approach one another as the separation increases, consistent with weak interactions.
Growth and relaxation to steady rotation
To compare directly with the material parameters used in Ref. [25 ] , the angular-speed calculations use
ϵ 21 = − 0.1097 , σ 21 ≃ − 0.5. (128) \epsilon_{21} = -0.1097,\qquad\sigma_{21} \simeq-0.5. \tag*{(128)} ϵ 21 = − 0.1097 , σ 21 ≃ − 0.5. (128) The field is E 0 / E c ( 0 ) = 2 E_0/E_c^{(0)} = 2 E 0 / E c ( 0 ) = 2 λ 0 = 3 \lambda_0 = 3 λ 0 = 3 3 \sqrt{3} 3 R = 4 , 6 , 8 R = 4, 6, 8 R = 4 , 6 , 8 Θ = 0 , π / 4 , π / 2 \Theta= 0, \pi/4, \pi/2 Θ = 0 , π /4 , π /2 X − X_- X − Y + Y_+ Y +
We compare the mean angular speed defined in Eq. (103) . The unexpanded equations are initialized about their interacting non-rotating state using the initial dipole perturbation measured from BEM. Time and angular velocity are expressed in the same Maxwell–Wagner units for all three descriptions.
Figure 11 shows the evolution of the mean angular speed from the perturbed non-rotating state. BEM and the unexpanded model exhibit the same rapid growth, overshoot, and damped approach to steady rotation. The horizontal dashed lines are the first-order regular asymptotic predictions in Eqs. (104) and (105) . They describe the final rotating state; the transient curves are obtained from BEM and the time-dependent unexpanded model. At the scale of Figure 11 , the nine histories have very similar growth, overshoot, and relaxation. Their steady speeds carry the smaller orientation-dependent differences: the BEM interaction-induced changes are at most about 1.5% of the matched isolated BEM speed. Table I resolves these steady differences.
Table I. Steady mean angular speeds at E 0 / E c ( 0 ) = 2 E_0/E_c^{(0)} = 2 E 0 / E c ( 0 ) = 2 (128) . The BEM and unexpanded values are averages after the transients have decayed; the final column gives the first-order regular asymptotic prediction. The last row gives the isolated references: matched BEM calculations give 1.73118, while the unexpanded model and regular asymptotics give 3 \sqrt{3} 3 τ M W − 1 \tau_{\mathrm{MW}}^{-1} τ MW − 1
Figure 11. Temporal evolution of the mean angular speed at E 0 / E c ( 0 ) = 2 E_0/E_{\mathrm{c}}^{(0)} = 2 E 0 / E c ( 0 ) = 2 (128) . Rows correspond to R = 4 , 6 , 8 R = 4, 6, 8 R = 4 , 6 , 8 Θ = 0 , π / 4 , π / 2 \Theta= 0, \pi/4, \pi/2 Θ = 0 , π /4 , π /2 X − X_{-} X − Y + Y_{+} Y + (103) . Time is scaled by τ M W \tau_{\mathrm{MW}} τ MW τ M W − 1 \tau_{\mathrm{MW}}^{-1} τ MW − 1
Effects of pair orientation and separation
The changes in steady rotation depend on the pair orientation. To measure the effect of the other sphere, we compare each pair with its corresponding isolated-sphere state. For BEM,
Δ Ω ‾ B E M = Ω ‾ p a i r , B E M − Ω ‾ i s o , B E M . (129) \overline{\Delta\Omega}_{\mathrm{BEM}} = \overline{\Omega}_{\mathrm{pair,BEM}} - \overline{\Omega}_{\mathrm{iso,BEM}}. \tag*{(129)} ΔΩ BEM = Ω pair , BEM − Ω iso , BEM . (129) The isolated calculation uses the same numerical settings as the pair calculation and gives Ω i s o , B E M ≃ 1.73118 \Omega_{\mathrm{iso,BEM}} \simeq1.73118 Ω iso , BEM ≃ 1.73118
3 ≃ 1.73205 \sqrt{3} \simeq1.73205 3 ≃ 1.73205 R = 8 R = 8 R = 8 Θ = π / 4 \Theta= \pi/4 Θ = π /4 3 \sqrt{3} 3
For this material and field strength, the interaction reduces the mean speed of the counter-rotating pair aligned with the field (Θ = 0 \Theta= 0 Θ = 0 Θ = π / 4 \Theta= \pi/4 Θ = π /4 Θ = π / 2 \Theta= \pi/2 Θ = π /2
The signs follow from the electric and hydrodynamic contributions in Eqs. (104) and (105) . For counter-rotation at Θ = 0 \Theta= 0 Θ = 0 Θ = π / 4 \Theta= \pi/4 Θ = π /4
At oblique pair angles, the two counter-rotating spheres can have different angular speeds. For example, at R = 4 R = 4 R = 4 Θ = π / 4 \Theta= \pi/4 Θ = π /4 e z e_z e z
As the separation increases, the interaction-induced changes become smaller, as expected from the R − 3 R^{-3} R − 3 Table I ). The additional mutual couplings in the unexpanded model produce only small changes to these steady speeds.
Conclusions
We have studied how electric and hydrodynamic interactions change the onset and rotation of two Quincke spheres at a prescribed pair configuration. At fixed field above onset, the first interaction correction to the steady angular speed is O ( R − 3 ) O(R^{-3}) O ( R − 3 ) λ 0 = O ( R − 3 ) \lambda_0 = O(R^{-3}) λ 0 = O ( R − 3 ) O ( R − 3 / 2 ) O(R^{-3/2}) O ( R − 3/2 )
The near-onset amplitude equations govern general transverse dipole perturbations. Their linearization identifies four collective modes, whose critical fields and associated single-mode rotating branches are determined together with their stability. Electric interaction gives a common leading shift of the modal thresholds, while hydrodynamic coupling splits them. Counter-rotation of type X − X_{-} X − 0 < Θ < Θ m 0 < \Theta< \Theta_{\mathrm{m}} 0 < Θ < Θ m Y + Y_{+} Y + Θ m < Θ ≤ π / 2 \Theta_{\mathrm{m}} < \Theta\le\pi/2 Θ m < Θ ≤ π /2 Θ m = cos − 1 ( 1 / 3 ) ≃ 54.7 ∘ \Theta_{\mathrm{m}} = \cos^{-1}(1/\sqrt{3}) \simeq54.7^\circ Θ m = cos − 1 ( 1/ 3 ) ≃ 54. 7 ∘ cos 2 Θ ∗ = [ 3 ( 1 + λ 0 ) ] − 1 \cos^2 \Theta_* = [3(1 + \lambda_0)]^{-1} cos 2 Θ ∗ = [ 3 ( 1 + λ 0 ) ] − 1 73.2 ∘ 73.2^\circ 73. 2 ∘ Θ m \Theta_{\mathrm{m}} Θ m δ ≪ λ 0 ≪ 1 \delta\ll\lambda_0 \ll1 δ ≪ λ 0 ≪ 1
Boundary element simulations provide onset intervals and follow the evolution of rotation from a small charge perturbation. The simulations and the unexpanded model show rapid growth, overshoot, and damped relaxation. For the tested material at twice the isolated critical field, the interaction slows a counter-rotating pair aligned with the field and speeds up the oblique counter-rotating and perpendicular co-rotating pairs. The two spheres in an oblique counter-rotating pair can also rotate at different angular speeds. These effects weaken as the separation grows. The unexpanded dipole–rotlet model provides a simple way to follow the rotation of both spheres by integrating six coupled equations for their dipole moments. It retains selected higher-order contributions and recovers the leading O ( R − 3 ) O(R^{-3}) O ( R − 3 ) Table I , its predictions of the steady mean angular speed agree with the BEM simulations to within 0.13%.
Allowing the particle centres to evolve with the computed translational velocities would couple rotation to changes in pair separation and orientation. The analysis can also be extended to larger collections of spheres, connecting with bulk suspension dynamics [59 ] , and to arrays of spatially constrained rotors motivated by recent micromotor experiments [56 ] . Other rigid shapes, including cylinders, would require their corresponding polarization and hydrodynamic responses. Incorporating wall-mediated interactions would permit a formal asymptotic treatment of pairs of Quincke rollers [41 , 42 ] . Further extensions include isolated deformable drops [74 ] , interacting drops [75 ] , and surfactant-covered drops [76 ] . For deformable drops, the solvability analysis must couple shape evolution to the interfacial flow, which depends on the viscosity ratio. For surfactant-covered drops, gradients in surfactant concentration also enter the tangential stress balance through Marangoni stresses.
Data availability
Supplementary Material is available. The data supporting this study are available from the corresponding author upon reasonable request.
Acknowledgments
This work was supported by the Leverhulme Trust through Research Project Grant No. RPG-2026-185.
References
Supplemental Material
Derivation of the dipole–rotlet equations
This section derives the unexpanded dipole–rotlet equations introduced in § 1.1 , Eqs. (7a)–(7c), of the main text. The surface-charge balance determines dipole evolution; the electric and hydrodynamic disturbances then couple the two spheres.
Surface-charge conservation
This subsection establishes the charge-conservation law underlying Eq. (7c) and the material charge balance in Eq. (124) of the main text. The corresponding node motion is given by Eq. (125) . We first use dimensional variables. A sphere of radius a a a x c ( t ) \mathbf{x}_c(t) x c ( t ) U \mathbf{U} U Ω \boldsymbol{\Omega} Ω s = x − x c \mathbf{s}=\mathbf{x}-\mathbf{x}_c s = x − x c r = ∣ s ∣ r=\lvert\mathbf{s}\rvert r = ∣ s ∣ n = s / a \mathbf{n}=\mathbf{s}/a n = s / a
V = U + Ω × s . (130) \mathbf{V}=\mathbf{U}+\boldsymbol{\Omega}\times\mathbf{s}. \tag*{(130)} V = U + Ω × s . (130) Let q s q_s q s
⟦ n ⋅ J ⟧ = n ⋅ ( J + − J − ) , J ± = σ 1 , 2 E ± . (131) \llbracket\mathbf{n}\cdot\mathbf{J}\rrbracket=\mathbf{n}\cdot\left(\mathbf{J}^{+}-\mathbf{J}^{-}\right),\qquad\mathbf{J}^{\pm}=\sigma_{1,2}\mathbf{E}^{\pm}. \tag*{(131)} [ [ n ⋅ J ] ] = n ⋅ ( J + − J − ) , J ± = σ 1 , 2 E ± . (131) The superscripts + + + − - −
D s q s D t + q s ∇ s ⋅ V + ⟦ n ⋅ J ⟧ = 0 , ∇ s = ( I − n n ) ⋅ ∇ . (132) \frac{D_s q_s}{Dt}+q_s\nabla_s\cdot\mathbf{V}+\llbracket\mathbf{n}\cdot\mathbf{J}\rrbracket=0,\qquad\nabla_s=\left(\mathbf{I}-\mathbf{n}\mathbf{n}\right)\cdot\nabla. \tag*{(132)} D t D s q s + q s ∇ s ⋅ V + [ [ n ⋅ J ] ] = 0 , ∇ s = ( I − nn ) ⋅ ∇. (132) Rigid motion preserves interfacial area, so ∇ s ⋅ V = 0 \nabla_s\cdot\mathbf{V}=0 ∇ s ⋅ V = 0
In coordinates translating with the centre but not rotating with the particle, define q ˉ s ( s , t ) = q s ( x c ( t ) + s , t ) \bar{q}_s(\mathbf{s},t)=q_s(\mathbf{x}_c(t)+\mathbf{s},t) q ˉ s ( s , t ) = q s ( x c ( t ) + s , t )
∂ t q ˉ s ∣ s = ∂ t q s ∣ x + U ⋅ ∇ q s . (133) \left.\partial_t\bar{q}_s\right|_{\mathbf{s}}=\left.\partial_tq_s\right|_{\mathbf{x}}+\mathbf{U}\cdot\nabla q_s. \tag*{(133)} ∂ t q ˉ s ∣ s = ∂ t q s ∣ x + U ⋅ ∇ q s . (133) Translation is therefore absorbed into the time derivative, and Eq. (132) becomes
∂ t q ˉ s + ( a Ω × n ) ⋅ ∇ s q ˉ s + ⟦ n ⋅ J ⟧ = 0. (134) \partial_t\bar{q}_s+(a\boldsymbol{\Omega}\times\mathbf{n})\cdot\nabla_s\bar{q}_s+\llbracket\mathbf{n}\cdot\mathbf{J}\rrbracket=0. \tag*{(134)} ∂ t q ˉ s + ( a Ω × n ) ⋅ ∇ s q ˉ s + [ [ n ⋅ J ] ] = 0. (134) The rotational surface velocity is divergence free, so the advection term can equivalently be written as a surface divergence. For the frozen pair geometry studied in the main text, the centres are not advanced and the retained surface motion is rotational.
Evolution of the electric dipole moment
We derive the dipole evolution law in Eq. (7c) of the main text by retaining the degree-one spherical harmonics of the charge and potential.
Potentials and boundary conditions
The boundary conditions used here are the same electric continuity and charge-jump conditions used in § 1.9 of the main text; the present calculation applies their dipolar projection. The bulk phases are charge free, giving ∇ 2 Φ ± = 0 \nabla^2\Phi^\pm=0 ∇ 2 Φ ± = 0 E ± = − ∇ Φ ± \mathbf{E}^\pm=-\nabla\Phi^\pm E ± = − ∇ Φ ± r = a r=a r = a
q s = n ⋅ ( ϵ 1 E + − ϵ 2 E − ) , (135) q_s=\mathbf{n}\cdot\left(\epsilon_1\mathbf{E}^+-\epsilon_2\mathbf{E}^-\right), \tag*{(135)} q s = n ⋅ ( ϵ 1 E + − ϵ 2 E − ) , (135) [ n ⋅ J ] = n ⋅ ( σ 1 E + − σ 2 E − ) . (136) [\mathbf{n}\cdot\mathbf{J}]=\mathbf{n}\cdot\left(\sigma_1\mathbf{E}^+-\sigma_2\mathbf{E}^-\right). \tag*{(136)} [ n ⋅ J ] = n ⋅ ( σ 1 E + − σ 2 E − ) . (136) Let E i ( t ) \mathbf{E}_i(t) E i ( t ) i i i E 0 e z E_0\mathbf{e}_z E 0 e z E i \mathbf{E}_i E i − r n ⋅ E i -r\mathbf{n}\cdot\mathbf{E}_i − r n ⋅ E i
Φ 1 + = − r n ⋅ E i + r − 2 n ⋅ P i , (137) \Phi_1^+=-r\mathbf{n}\cdot\mathbf{E}_i+r^{-2}\mathbf{n}\cdot\mathbf{P}_i, \tag*{(137)} Φ 1 + = − r n ⋅ E i + r − 2 n ⋅ P i , (137) Φ 1 − = − r n ⋅ E i + r a 3 n ⋅ P i . (138) \Phi_1^-=-r\mathbf{n}\cdot\mathbf{E}_i+\frac{r}{a^3}\mathbf{n}\cdot\mathbf{P}_i. \tag*{(138)} Φ 1 − = − r n ⋅ E i + a 3 r n ⋅ P i . (138) Here P i \mathbf{P}_i P i a 3 E 0 a^3E_0 a 3 E 0 § 1.1 , gives the dimensionless P i \mathbf{P}_i P i
Charge density and normal current
This step expresses the two terms in the charge balance using the material contrasts defined in Eq. (1) . The normal derivatives at the interface are
∂ r Φ 1 + ∣ r = a = − n ⋅ E i − 2 a 3 n ⋅ P i , (139) \left.\partial_r\Phi_1^+\right|_{r=a}=-\mathbf{n}\cdot\mathbf{E}_i-\frac{2}{a^3}\mathbf{n}\cdot\mathbf{P}_i, \tag*{(139)} ∂ r Φ 1 + r = a = − n ⋅ E i − a 3 2 n ⋅ P i , (139) ∂ r Φ 1 − ∣ r = a = − n ⋅ E i + 1 a 3 n ⋅ P i . (140) \left.\partial_r\Phi_1^-\right|_{r=a}=-\mathbf{n}\cdot\mathbf{E}_i+\frac{1}{a^3}\mathbf{n}\cdot\mathbf{P}_i. \tag*{(140)} ∂ r Φ 1 − r = a = − n ⋅ E i + a 3 1 n ⋅ P i . (140) Substitution into Eqs. (135) and (S7) gives
q s , 1 = ϵ 2 + 2 ϵ 1 a 3 n ⋅ ( P i − a 3 ϵ 21 E i ) , (141) q_{s,1}=\frac{\epsilon_2+2\epsilon_1}{a^3}\mathbf{n}\cdot\left(\mathbf{P}_i-a^3\epsilon_{21}\mathbf{E}_i\right), \tag*{(141)} q s , 1 = a 3 ϵ 2 + 2 ϵ 1 n ⋅ ( P i − a 3 ϵ 21 E i ) , (141) ( [ n ⋅ J ] ) 1 = σ 2 + 2 σ 1 a 3 n ⋅ ( P i − a 3 σ 21 E i ) , (142) \left([\mathbf{n}\cdot\mathbf{J}]\right)_1=\frac{\sigma_2+2\sigma_1}{a^3}\mathbf{n}\cdot\left(\mathbf{P}_i-a^3\sigma_{21}\mathbf{E}_i\right), \tag*{(142)} ( [ n ⋅ J ] ) 1 = a 3 σ 2 + 2 σ 1 n ⋅ ( P i − a 3 σ 21 E i ) , (142) where
ϵ 21 = ϵ 2 − ϵ 1 ϵ 2 + 2 ϵ 1 , σ 21 = σ 2 − σ 1 σ 2 + 2 σ 1 , τ M W = ϵ 2 + 2 ϵ 1 σ 2 + 2 σ 1 . (143) \epsilon_{21}=\frac{\epsilon_2-\epsilon_1}{\epsilon_2+2\epsilon_1},\qquad
\sigma_{21}=\frac{\sigma_2-\sigma_1}{\sigma_2+2\sigma_1},\qquad
\tau_{\mathrm{MW}}=\frac{\epsilon_2+2\epsilon_1}{\sigma_2+2\sigma_1}. \tag*{(143)} ϵ 21 = ϵ 2 + 2 ϵ 1 ϵ 2 − ϵ 1 , σ 21 = σ 2 + 2 σ 1 σ 2 − σ 1 , τ MW = σ 2 + 2 σ 1 ϵ 2 + 2 ϵ 1 . (143) The permittivity-controlled polarization enters the free-charge density, whereas the conductivity-controlled polarization enters the Ohmic current. This distinction determines the two different field terms in the dipole evolution equation.
Rigid rotation of the dipolar charge
The following identity supplies the rotational convection term in Eq. (7c) of the main text. For a vector B ( t ) \mathbf{B}(t) B ( t )
∇ s ( n ⋅ B ) = 1 a ( I − n n ) ⋅ B . (144) \nabla_s(\mathbf{n}\cdot\mathbf{B})=\frac{1}{a}(\mathbf{I}-\mathbf{n}\mathbf{n})\cdot\mathbf{B}. \tag*{(144)} ∇ s ( n ⋅ B ) = a 1 ( I − nn ) ⋅ B . (144) Consequently,
( a Ω i × n ) ⋅ ∇ s ( n ⋅ B ) = − n ⋅ ( Ω i × B ) . (145) (a\boldsymbol{\Omega}_i\times\mathbf{n})\cdot\nabla_s(\mathbf{n}\cdot\mathbf{B})=-\mathbf{n}\cdot(\boldsymbol{\Omega}_i\times\mathbf{B}). \tag*{(145)} ( a Ω i × n ) ⋅ ∇ s ( n ⋅ B ) = − n ⋅ ( Ω i × B ) . (145) Rigid rotation preserves the spherical-harmonic degree, so no higher harmonic is generated by advecting q s , 1 q_{s,1} q s , 1 B = P i − a 3 ϵ 21 E i \mathbf{B}=\mathbf{P}_i-a^3\epsilon_{21}\mathbf{E}_i B = P i − a 3 ϵ 21 E i (145) evaluates the rotational term in Eq. (134) .
Dipole moment equation in a time-dependent incident field
We now assemble the dimensional dipole equation and recover Eq. (7c), including the local-field derivative. At fixed s \boldsymbol{s} s (141) gives
∂ t q s , 1 = ϵ 2 + 2 ϵ 1 a 3 n ⋅ ( P ˙ i − a 3 ϵ 21 E ˙ i ) . (146) \partial_t q_{s,1}=\frac{\epsilon_2+2\epsilon_1}{a^3}\,\boldsymbol{n}\cdot\left(\dot{\boldsymbol{P}}_i-a^3\epsilon_{21}\dot{\boldsymbol{E}}_i\right). \tag*{(146)} ∂ t q s , 1 = a 3 ϵ 2 + 2 ϵ 1 n ⋅ ( P ˙ i − a 3 ϵ 21 E ˙ i ) . (146) Combining this result with Eqs. (142) , (S16), and (S5) gives a scalar identity valid for every n \boldsymbol{n} n ∫ S 2 n n d Ω = ( 4 π / 3 ) I \int_{S^2}\boldsymbol{n}\boldsymbol{n}\,\mathrm{d}\Omega=(4\pi/3)\boldsymbol{I} ∫ S 2 n n d Ω = ( 4 π /3 ) I
P ˙ i = Ω i × ( P i − a 3 ϵ 21 E i ) − 1 τ M W ( P i − a 3 σ 21 E i ) + a 3 ϵ 21 E ˙ i . (147) \dot{\boldsymbol{P}}_i=\boldsymbol{\Omega}_i\times\left(\boldsymbol{P}_i-a^3\epsilon_{21}\boldsymbol{E}_i\right)-\frac{1}{\tau_{\mathrm{MW}}}\left(\boldsymbol{P}_i-a^3\sigma_{21}\boldsymbol{E}_i\right)+a^3\epsilon_{21}\dot{\boldsymbol{E}}_i. \tag*{(147)} P ˙ i = Ω i × ( P i − a 3 ϵ 21 E i ) − τ MW 1 ( P i − a 3 σ 21 E i ) + a 3 ϵ 21 E ˙ i . (147) The last term arises from differentiating the incident-field contribution to the free charge. The other sphere’s evolving dipole makes this local field time dependent.
Scaling as in the main text removes a a a τ M W \tau_{\mathrm{MW}} τ MW
P ˙ i = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) + ϵ 21 E ˙ i . (148) \dot{\boldsymbol{P}}_i=\boldsymbol{\Omega}_i\times\left(\boldsymbol{P}_i-\epsilon_{21}\boldsymbol{E}_i\right)-\left(\boldsymbol{P}_i-\sigma_{21}\boldsymbol{E}_i\right)+\epsilon_{21}\dot{\boldsymbol{E}}_i. \tag*{(148)} P ˙ i = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) + ϵ 21 E ˙ i . (148) The isolated critical Mason number is M a c \mathrm{Ma}_{\mathrm{c}} Ma c (3) . In dimensional variables this corresponds to
[ E c ( 0 ) ] 2 = 2 μ τ M W ϵ 1 M a c . (149) \left[E_{\mathrm{c}}^{(0)}\right]^2=\frac{2\mu}{\tau_{\mathrm{MW}}\epsilon_1\mathrm{Ma}_{\mathrm{c}}}. \tag*{(149)} [ E c ( 0 ) ] 2 = τ MW ϵ 1 Ma c 2 μ . (149) The required material condition ϵ 21 > σ 21 \epsilon_{21}>\sigma_{21} ϵ 21 > σ 21 ϵ 2 σ 1 > ϵ 1 σ 2 \epsilon_2\sigma_1>\epsilon_1\sigma_2 ϵ 2 σ 1 > ϵ 1 σ 2
Interaction fields and coupled dipoles
This subsection derives the disturbance fields in Eqs. (7a)–(7c) of the main text.
Electric dipole interaction
The electric disturbance gives Eq. (7a). In the dimensionless variables of the main text, a dipole P j \boldsymbol{P}_j P j x c , j \boldsymbol{x}_{c,j} x c , j
Φ j d ( x ) = s j ⋅ P j s j 3 , s j = x − x c , j , s j = ∣ s j ∣ . (150) \Phi_j^{d}(\boldsymbol{x})=\frac{\boldsymbol{s}_j\cdot\boldsymbol{P}_j}{s_j^3},\qquad\boldsymbol{s}_j=\boldsymbol{x}-\boldsymbol{x}_{c,j},\qquad s_j=\lvert\boldsymbol{s}_j\rvert. \tag*{(150)} Φ j d ( x ) = s j 3 s j ⋅ P j , s j = x − x c , j , s j = ∣ s j ∣ . (150) Its negative gradient is
E j d = − s j − 3 ( I − 3 s ^ j s ^ j ) ⋅ P j . (151) \boldsymbol{E}_j^{d}=-s_j^{-3}\left(\boldsymbol{I}-3\hat{\boldsymbol{s}}_j\hat{\boldsymbol{s}}_j\right)\cdot\boldsymbol{P}_j. \tag*{(151)} E j d = − s j − 3 ( I − 3 s ^ j s ^ j ) ⋅ P j . (151) At the other centre, the local field is therefore
E i = e z − δ Π ⋅ P j , Π = I − 3 R ^ R ^ , δ = R − 3 . (152) \boldsymbol{E}_i=\boldsymbol{e}_z-\delta\boldsymbol{\Pi}\cdot\boldsymbol{P}_j,\qquad\boldsymbol{\Pi}=\boldsymbol{I}-3\hat{\boldsymbol{R}}\hat{\boldsymbol{R}},\qquad\delta=R^{-3}. \tag*{(152)} E i = e z − δ Π ⋅ P j , Π = I − 3 R ^ R ^ , δ = R − 3 . (152) The dyadic tensor is unchanged on reversing the pair axis, so the same expression applies to both particles. It has eigenvalues − 2 -2 − 2 1 1 1
Rotlet interaction and torque balance
The rotational response to the fluid disturbance gives Eq. (7b). The leading disturbance of a rotating unit sphere is
u j d = Ω j × s j s j 3 . (153) \boldsymbol{u}_j^{d}=\frac{\boldsymbol{\Omega}_j\times\boldsymbol{s}_j}{s_j^3}. \tag*{(153)} u j d = s j 3 Ω j × s j . (153) For a constant vector Ω j \boldsymbol{\Omega}_j Ω j s j ≠ 0 s_j \ne0 s j = 0
∇ × u j d = − s j − 3 ( I − 3 s ^ j s ^ j ) ⋅ Ω j . (154) \nabla\times\mathbf{u}_j^d=-s_j^{-3}\left(\mathbf{I}-3\hat{\mathbf{s}}_j\hat{\mathbf{s}}_j\right)\cdot\boldsymbol{\Omega}_j. \tag*{(154)} ∇ × u j d = − s j − 3 ( I − 3 s ^ j s ^ j ) ⋅ Ω j . (154) The electric torque on a sphere in a locally uniform field is 4 π ϵ 1 P i × E i 4\pi\epsilon_1\mathbf{P}_i \times\mathbf{E}_i 4 π ϵ 1 P i × E i − 8 π μ a 3 [ Ω i − 1 2 ∇ × u e ] -8\pi\mu a^3[\boldsymbol{\Omega}_i-\frac{1}{2}\nabla\times\mathbf{u}_e] − 8 π μ a 3 [ Ω i − 2 1 ∇ × u e ] u e \mathbf{u}_e u e i i i u j d \mathbf{u}_j^d u j d (153) , with j ≠ i j\ne i j = i
Ω i + δ 2 Π ⋅ Ω j = M a − 1 P i × E i . (155) \boldsymbol{\Omega}_i+\frac{\delta}{2}\mathbf{\Pi}\cdot\boldsymbol{\Omega}_j=\mathrm{Ma}^{-1}\mathbf{P}_i\times\mathbf{E}_i. \tag*{(155)} Ω i + 2 δ Π ⋅ Ω j = Ma − 1 P i × E i . (155) The source strength in Eq. (153) is the isolated-sphere angular velocity at leading interaction order. Retaining the coupled angular velocities in Eq. (155) defines the unexpanded dipole–rotlet model used here.
Coupled dipole derivatives
The local-field derivative couples the two dipole rates in Eq. (7c). For fixed centres, δ \delta δ Π \mathbf{\Pi} Π (152) gives
E ˙ i = − δ Π ⋅ P ˙ j . (156) \dot{\mathbf{E}}_i=-\delta\mathbf{\Pi}\cdot\dot{\mathbf{P}}_j. \tag*{(156)} E ˙ i = − δ Π ⋅ P ˙ j . (156) The nondimensional dipole equation (S19) can therefore be written as
P ˙ i + ϵ 21 δ Π ⋅ P ˙ j = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) . (157) \dot{\mathbf{P}}_i+\epsilon_{21}\delta\mathbf{\Pi}\cdot\dot{\mathbf{P}}_j=\boldsymbol{\Omega}_i\times\left(\mathbf{P}_i-\epsilon_{21}\mathbf{E}_i\right)-\left(\mathbf{P}_i-\sigma_{21}\mathbf{E}_i\right). \tag*{(157)} P ˙ i + ϵ 21 δ Π ⋅ P ˙ j = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) . (157) Equations (S23), (S26), and (S28) are the unexpanded dipole–rotlet model used in Eq. (7) of the main text.
To implement the model, let f i = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i ) f_i=\boldsymbol{\Omega}_i\times(\mathbf{P}_i-\epsilon_{21}\mathbf{E}_i)-(\mathbf{P}_i-\sigma_{21}\mathbf{E}_i) f i = Ω i × ( P i − ϵ 21 E i ) − ( P i − σ 21 E i )
\begin{pmatrix}
\mathbf{I} & \delta\mathbf{\Pi}/2 \\
\delta\mathbf{\Pi}/2 & \mathbf{I}
\end{pmatrix}
\begin{pmatrix}
\boldsymbol{\Omega}_1 \\
\boldsymbol{\Omega}_2
\end{pmatrix}
=\mathrm{Ma}^{-1}
\begin{pmatrix}
\mathbf{P}_1\times\mathbf{E}_1 \\
\mathbf{P}_2\times\mathbf{E}_2
\end{pmatrix}. \tag*{(158)}\end{pmatrix}
\begin{pmatrix}
\mathbf{I} & \epsilon_{21}\delta\mathbf{\Pi} \\
\epsilon_{21}\delta\mathbf{\Pi} & \mathbf{I}
\end{pmatrix}
\begin{pmatrix}
\dot{\mathbf{P}}_1 \\
\dot{\mathbf{P}}_2
\end{pmatrix}
=
\begin{pmatrix}
\mathbf{f}_1 \\
\mathbf{f}_2
\end{pmatrix}. \tag*{(159)}\end{pmatrix}
The torque matrix has eigenvalues 1 ± δ k / 2 1\pm\delta k/2 1 ± δ k /2 1 ± ϵ 21 δ k 1\pm\epsilon_{21}\delta k 1 ± ϵ 21 δ k k ∈ { − 2 , 1 , 1 } k\in\{-2,1,1\} k ∈ { − 2 , 1 , 1 } − 1 / 2 < ϵ 21 < 1 -1/2<\epsilon_{21}<1 − 1/2 < ϵ 21 < 1 δ < 1 / 8 \delta<1/8 δ < 1/8
Phase equilibria and stability away from onset
This section supplements Sec. IV D and Figure 2 of the main text. The counter-rotating and co-rotating stability conditions are derived there from Eq. (38) . Here we examine the additional branch with ϕ 2 = − ϕ 1 \phi_2=-\phi_1 ϕ 2 = − ϕ 1 δ / λ 0 ≪ 1 \delta/\lambda_0\ll1 δ / λ 0 ≪ 1 Π x x = 3 cos 2 Θ − 2 \Pi_{xx}=3\cos^2\Theta-2 Π xx = 3 cos 2 Θ − 2 Π z x = − 3 sin Θ cos Θ \Pi_{zx}=-3\sin\Theta\cos\Theta Π z x = − 3 sin Θ cos Θ
To examine the degeneracy, define the half-sum and half-difference ϕ s = ( ϕ 1 + ϕ 2 ) / 2 \phi_s=(\phi_1+\phi_2)/2 ϕ s = ( ϕ 1 + ϕ 2 ) /2 ϕ a = ( ϕ 1 − ϕ 2 ) / 2 \phi_a=(\phi_1-\phi_2)/2 ϕ a = ( ϕ 1 − ϕ 2 ) /2 (38) becomes
d ϕ s d τ = sin ϕ s 2 [ ( 1 − Π x x ) cos ϕ s − λ 0 Π z x cos ϕ a ] , (160) \frac{\mathrm{d}\phi_s}{\mathrm{d}\tau}=\frac{\sin\phi_s}{2}\left[(1-\Pi_{xx})\cos\phi_s-\sqrt{\lambda_0}\Pi_{zx}\cos\phi_a\right], \tag*{(160)} d τ d ϕ s = 2 sin ϕ s [ ( 1 − Π xx ) cos ϕ s − λ 0 Π z x cos ϕ a ] , (160) d ϕ a d τ = sin ϕ a 2 [ ( 1 + Π x x ) cos ϕ a + λ 0 Π z x cos ϕ s ] . (161) \frac{\mathrm{d}\phi_a}{\mathrm{d}\tau}=\frac{\sin\phi_a}{2}\left[(1+\Pi_{xx})\cos\phi_a+\sqrt{\lambda_0}\Pi_{zx}\cos\phi_s\right]. \tag*{(161)} d τ d ϕ a = 2 sin ϕ a [ ( 1 + Π xx ) cos ϕ a + λ 0 Π z x cos ϕ s ] . (161) The mirror branch ϕ 2 = − ϕ 1 \phi_2=-\phi_1 ϕ 2 = − ϕ 1 ϕ s = 0 \phi_s=0 ϕ s = 0 cos ϕ a = − λ 0 Π z x / ( 1 + Π x x ) \cos\phi_a=-\sqrt{\lambda_0}\Pi_{zx}/(1+\Pi_{xx}) cos ϕ a = − λ 0 Π z x / ( 1 + Π xx ) ( 1 − Π x x 2 + λ 0 Π z x 2 ) / [ 2 ( 1 + Π x x ) ] (1-\Pi_{xx}^2+\lambda_0\Pi_{zx}^2)/[2(1+\Pi_{xx})] ( 1 − Π xx 2 + λ 0 Π z x 2 ) / [ 2 ( 1 + Π xx )] − ( 1 + Π x x ) sin 2 ϕ a / 2 -(1+\Pi_{xx})\sin^2\phi_a/2 − ( 1 + Π xx ) sin 2 ϕ a /2 Θ ∗ \Theta_\ast Θ ∗ (160) and (S32) become dependent:
( 1 − Π x x ) cos ϕ s = λ 0 Π z x cos ϕ a , 1 − Π x x 2 + λ 0 Π z x 2 = 0. (162) \begin{aligned}
(1-\Pi_{xx})\cos\phi_s&=\sqrt{\lambda_0}\Pi_{zx}\cos\phi_a,\\
1-\Pi_{xx}^2+\lambda_0\Pi_{zx}^2&=0. \tag*{(162)}
\end{aligned} ( 1 − Π xx ) cos ϕ s 1 − Π xx 2 + λ 0 Π z x 2 = λ 0 Π z x cos ϕ a , = 0. (162) They give a continuous family connecting the counter-rotating and co-rotating states. As λ 0 → 0 + \lambda_0 \to0^+ λ 0 → 0 + ϕ s = π / 2 \phi_s = \pi/2 ϕ s = π /2 X − X_- X − Y + Y_+ Y +
Stability of the near-onset rotating branches
This section supplements Sec. VII B and Figure 7 (b) of the main text. The order-by-order solvability calculation is given in Sec. VI; here we give the complete linear stability calculation for the four single-mode branches. The X − X_- X − (88) of the main text. Let z i = ( X i , Y i ) T \mathbf{z}_i=(X_i,Y_i)^{\mathsf{T}} z i = ( X i , Y i ) T g = ρ − 2 σ 21 ( 1 − 3 cos 2 Θ ) g=\rho-2\sigma_{21}(1-3\cos^2\Theta) g = ρ − 2 σ 21 ( 1 − 3 cos 2 Θ ) (76) and (77) of the main text can be written as
d z i d T 2 = g z i + K z j − M a c − 2 ∣ z i ∣ 2 z i , K = diag ( − 1 / 2 , h ) . (163) \frac{\mathrm{d}\mathbf{z}_i}{\mathrm{d}T_2}=g\mathbf{z}_i+K\mathbf{z}_j-Ma_c^{-2}\lvert\mathbf{z}_i\rvert^2\mathbf{z}_i,\qquad K=\operatorname{diag}(-1/2,h). \tag*{(163)} d T 2 d z i = g z i + K z j − M a c − 2 ∣ z i ∣ 2 z i , K = diag ( − 1/2 , h ) . (163) Here h = ( 2 − 3 cos 2 Θ ) / 2 h=(2-3\cos^2\Theta)/2 h = ( 2 − 3 cos 2 Θ ) /2 (78) of the main text. At any equilibrium, the Jacobian in the ordering ( X 1 , Y 1 , X 2 , Y 2 ) (X_1,Y_1,X_2,Y_2) ( X 1 , Y 1 , X 2 , Y 2 )
J ∗ = ( A 1 K K A 2 ) , A i = ( g − M a c − 2 ∣ z i ∣ 2 ) I 2 − 2 M a c − 2 z i z i T . (164) J^*=\begin{pmatrix} A_1 & K \\ K & A_2 \end{pmatrix},\qquad A_i=(g-Ma_c^{-2}\lvert\mathbf{z}_i\rvert^2)I_2-2Ma_c^{-2}\mathbf{z}_i\mathbf{z}_i^{\mathsf{T}}. \tag*{(164)} J ∗ = ( A 1 K K A 2 ) , A i = ( g − M a c − 2 ∣ z i ∣ 2 ) I 2 − 2 M a c − 2 z i z i T . (164) For a state with only an X X X Y Y Y (80) , now for infinitesimal perturbations of a rotating branch:
Δ X ± = Δ X 1 ± Δ X 2 2 , Δ Y ± = Δ Y 1 ± Δ Y 2 2 . (165) \Delta X_{\pm}=\frac{\Delta X_1\pm\Delta X_2}{\sqrt{2}},\qquad\Delta Y_{\pm}=\frac{\Delta Y_1\pm\Delta Y_2}{\sqrt{2}}. \tag*{(165)} Δ X ± = 2 Δ X 1 ± Δ X 2 , Δ Y ± = 2 Δ Y 1 ± Δ Y 2 . (165) The columns of Table S1 give the growth rates of equal (+) or opposite (−) perturbations on the two spheres in the X X X Y Y Y T 2 T_2 T 2
For a unique largest mode shift, the corresponding nonzero branch has negative magnitude and direction growth rates. When the line of centres is parallel to the field, the largest shifts m X − m_{X_-} m X − m Y − m_{Y_-} m Y −
At Θ = Θ m \Theta=\Theta_m Θ = Θ m cos 2 Θ = 1 / 3 \cos^2\Theta=1/3 cos 2 Θ = 1/3 X + = Y − = 0 X_+=Y_-=0 X + = Y − = 0
X − 2 + Y + 2 = 2 M a c 2 ( ρ + 1 2 ) , ρ > − 1 2 . (166) X_-^2+Y_+^2=2Ma_c^2\left(\rho+\frac{1}{2}\right),\qquad\rho>-\frac{1}{2}. \tag*{(166)} X − 2 + Y + 2 = 2 M a c 2 ( ρ + 2 1 ) , ρ > − 2 1 . (166) Thus a continuous family with both X − X_- X − Y + Y_+ Y +
Critical fields and steady angular velocities
This section supplies the supporting calculations for the unexpanded stationary thresholds in Sec. VIII and the regular asymptotic angular speeds in Sec. IX of the main text.
Stationary critical fields of the unexpanded model
We derive and justify the stationary onset fields in Sec. VIII B, Eqs. (96) and (97) , before evaluating their crossing.
Linearized equations for equal and opposite perturbations
The linearized dipole equations below supplement the neutral angular-velocity equation (95) of the main text. We retain the coupled dipole–rotlet equations without expanding in δ \delta δ (98) and (99) of the main text. In Eq. (89) of the main text, the switches s E s_E s E s H s_H s H
E n r = ( I + s E σ 21 δ Π ) − 1 ⋅ e z , P n r = σ 21 E n r . (167) \mathbf{E}_{\mathrm{nr}} = \left(I+s_E\sigma_{21}\delta\Pi\right)^{-1}\cdot\mathbf{e}_z,\qquad\mathbf{P}_{\mathrm{nr}}=\sigma_{21}\mathbf{E}_{\mathrm{nr}}. \tag*{(167)} E nr = ( I + s E σ 21 δ Π ) − 1 ⋅ e z , P nr = σ 21 E nr . (167) Write P i = P n r + p i \mathbf{P}_i=\mathbf{P}_{\mathrm{nr}}+\mathbf{p}_i P i = P nr + p i E i = E n r + e i \mathbf{E}_i=\mathbf{E}_{\mathrm{nr}}+\mathbf{e}_i E i = E nr + e i Ω i = ω i \boldsymbol{\Omega}_i=\boldsymbol{\omega}_i Ω i = ω i e i = − s E δ Π ⋅ p j \mathbf{e}_i=-s_E\delta\Pi\cdot\mathbf{p}_j e i = − s E δ Π ⋅ p j
ω i + s H δ 2 Π ⋅ ω j = M a − 1 ( p i × E n r + P n r × e i ) , (168) \boldsymbol{\omega}_i+s_H\frac{\delta}{2}\Pi\cdot\boldsymbol{\omega}_j=\mathrm{Ma}^{-1}\left(\mathbf{p}_i\times\mathbf{E}_{\mathrm{nr}}+\mathbf{P}_{\mathrm{nr}}\times\mathbf{e}_i\right), \tag*{(168)} ω i + s H 2 δ Π ⋅ ω j = Ma − 1 ( p i × E nr + P nr × e i ) , (168) p ˙ i + s E ϵ 21 δ Π ⋅ p ˙ j = M a c C ( E n r ) ⋅ ω i − p i + σ 21 e i . (169) \dot{\mathbf{p}}_i+s_E\epsilon_{21}\delta\Pi\cdot\dot{\mathbf{p}}_j=\mathrm{Ma}_c\mathcal{C}(\mathbf{E}_{\mathrm{nr}})\cdot\boldsymbol{\omega}_i-\mathbf{p}_i+\sigma_{21}\mathbf{e}_i. \tag*{(169)} p ˙ i + s E ϵ 21 δ Π ⋅ p ˙ j = Ma c C ( E nr ) ⋅ ω i − p i + σ 21 e i . (169) Here C ( v ) w = v × w \mathcal{C}(\mathbf{v})\mathbf{w}=\mathbf{v}\times\mathbf{w} C ( v ) w = v × w p 2 = q p 1 \mathbf{p}_2=q\mathbf{p}_1 p 2 = q p 1 ω 2 = q ω 1 \boldsymbol{\omega}_2=q\boldsymbol{\omega}_1 ω 2 = q ω 1 q = ± 1 q=\pm1 q = ± 1 p q = p 1 \mathbf{p}_q=\mathbf{p}_1 p q = p 1 e q = e 1 \mathbf{e}_q=\mathbf{e}_1 e q = e 1 ω q = ω 1 \boldsymbol{\omega}_q=\boldsymbol{\omega}_1 ω q = ω 1
B q = I + q s E σ 21 δ Π , M q = I + q s E ϵ 21 δ Π , H q = I + q s H δ 2 Π . (170) B_q=I+qs_E\sigma_{21}\delta\Pi,\qquad M_q=I+qs_E\epsilon_{21}\delta\Pi,\qquad H_q=I+qs_H\frac{\delta}{2}\Pi. \tag*{(170)} B q = I + q s E σ 21 δ Π , M q = I + q s E ϵ 21 δ Π , H q = I + q s H 2 δ Π. (170) Eliminating angular velocity gives
M q p ˙ q = ( ( E 0 E c ( 0 ) ) 2 L q − I ) B q p q , M a c M a = ( E 0 E c ( 0 ) ) 2 , (171) M_q\dot{\mathbf{p}}_q=\left(\left(\frac{E_0}{E_c^{(0)}}\right)^2L_q-I\right)B_q\mathbf{p}_q,\qquad\frac{\mathrm{Ma}_c}{\mathrm{Ma}}=\left(\frac{E_0}{E_c^{(0)}}\right)^2, \tag*{(171)} M q p ˙ q = ( ( E c ( 0 ) E 0 ) 2 L q − I ) B q p q , Ma Ma c = ( E c ( 0 ) E 0 ) 2 , (171) L q = − C ( E n r ) H q − 1 C ( E n r ) . (172) L_q=-\mathcal{C}(\mathbf{E}_{\mathrm{nr}})H_q^{-1}\mathcal{C}(\mathbf{E}_{\mathrm{nr}}). \tag*{(172)} L q = − C ( E nr ) H q − 1 C ( E nr ) . (172) The two exact three-dimensional Jacobians are therefore
J q = M q − 1 ( ( E 0 E c ( 0 ) ) 2 L q − I ) B q . (173) J_q=M_q^{-1}\left(\left(\frac{E_0}{E_c^{(0)}}\right)^2L_q-I\right)B_q. \tag*{(173)} J q = M q − 1 ( ( E c ( 0 ) E 0 ) 2 L q − I ) B q . (173) In the particle ordering ( p 1 , p 2 ) (\mathbf{p}_1,\mathbf{p}_2) ( p 1 , p 2 )
J 6 = 1 2 ( J + + J − J + − J − J + − J − J + + J − ) . (174) J_6=\frac{1}{2}\begin{pmatrix}
J_++J_- & J_+-J_-\\
J_+-J_- & J_++J_-
\end{pmatrix}. \tag*{(174)} J 6 = 2 1 ( J + + J − J + − J − J + − J − J + + J − ) . (174) Reality of the growth rates and stationary onset
This subsection establishes the reality of the growth rates underlying the stationary-onset calculation in Sec. VIII B. For sufficiently small δ \delta δ B q B_q B q M q M_q M q H q H_q H q Π \Pi Π u q = B q p q = p q − σ 21 e q \mathbf{u}_q=B_q\mathbf{p}_q=\mathbf{p}_q-\sigma_{21}\mathbf{e}_q u q = B q p q = p q − σ 21 e q e q = − q s E δ Π ⋅ p q \mathbf{e}_q=-qs_E\delta\Pi\cdot\mathbf{p}_q e q = − q s E δ Π ⋅ p q u q \mathbf{u}_q u q
D q u ˙ q = ( ( E 0 E c ( 0 ) ) 2 L q − I ) u q , D q = M q B q − 1 . (175) D_q\dot{\mathbf{u}}_q=\left(\left(\frac{E_0}{E_c^{(0)}}\right)^2L_q-I\right)\mathbf{u}_q,\qquad D_q=M_qB_q^{-1}. \tag*{(175)} D q u ˙ q = ( ( E c ( 0 ) E 0 ) 2 L q − I ) u q , D q = M q B q − 1 . (175) The mass matrix D q D_q D q D q − 1 / 2 ( ( E 0 E c ( 0 ) ) 2 L q − I ) D q − 1 / 2 D_q^{-1/2}\left(\left(\frac{E_0}{E_c^{(0)}}\right)^2L_q-I\right)D_q^{-1/2} D q − 1/2 ( ( E c ( 0 ) E 0 ) 2 L q − I ) D q − 1/2 L q L_q L q C T = − C C^{\mathsf{T}}=-C C T = − C H q H_q H q E n r \mathbf{E}_{\mathrm{nr}} E nr E 0 / E c ( 0 ) = 0 E_0/E_c^{(0)}=0 E 0 / E c ( 0 ) = 0 ( E 0 / E c ( 0 ) ) 2 \left(E_0/E_c^{(0)}\right)^2 ( E 0 / E c ( 0 ) ) 2 L q L_q L q
Since B q B_q B q
u q = ( E 0 E c ( 0 ) ) 2 L q u q , u q = B q p q . (176) \mathbf{u}_q=\left(\frac{E_0}{E_c^{(0)}}\right)^2L_q\mathbf{u}_q,\qquad\mathbf{u}_q=B_q\mathbf{p}_q. \tag*{(176)} u q = ( E c ( 0 ) E 0 ) 2 L q u q , u q = B q p q . (176) Neither M q M_q M q B q B_q B q ϵ 21 \epsilon_{21} ϵ 21
Calculation of the four stationary critical fields
This calculation evaluates the two mode directions for each exchange symmetry, yielding Eqs. (96) and (97) of the main text. Take R ^ \hat{\mathbf{R}} R ^ E n r \mathbf{E}_{\mathrm{nr}} E nr x x x z z z
F = ∣ E n r ∣ 2 , G = ( R ^ ⋅ E n r ) 2 . (177) F = \lvert\mathbf{E}_{\mathrm{nr}}\rvert^2,\qquad G = \left(\hat{\mathbf{R}}\cdot\mathbf{E}_{\mathrm{nr}}\right)^2. \tag*{(177)} F = ∣ E nr ∣ 2 , G = ( R ^ ⋅ E nr ) 2 . (177) The X q X_q X q L q \mathbf{L}_q L q E n r \mathbf{E}_{\mathrm{nr}} E nr E n r \mathbf{E}_{\mathrm{nr}} E nr e y \mathbf{e}_y e y H q \mathbf{H}_q H q 1 + q s H δ / 2 1+q s_{\mathrm{H}}\delta/2 1 + q s H δ /2 Y q Y_q Y q e y \mathbf{e}_y e y
ℓ X q = F 1 + q s H δ / 2 , ℓ Y q = ( 1 + q s H δ / 2 ) F − ( 3 q s H δ / 2 ) G ( 1 + q s H δ / 2 ) ( 1 − q s H δ ) . (178) \ell_{X_q} = \frac{F}{1+q s_{\mathrm{H}}\delta/2},\qquad
\ell_{Y_q} = \frac{(1+q s_{\mathrm{H}}\delta/2)F-(3q s_{\mathrm{H}}\delta/2)G}{(1+q s_{\mathrm{H}}\delta/2)(1-q s_{\mathrm{H}}\delta)}. \tag*{(178)} ℓ X q = 1 + q s H δ /2 F , ℓ Y q = ( 1 + q s H δ /2 ) ( 1 − q s H δ ) ( 1 + q s H δ /2 ) F − ( 3 q s H δ /2 ) G . (178) Taking their reciprocals gives Eqs. (96) and (97) of the main text, which are exact within the unexpanded model. The in-plane dipole eigenvector is obtained by applying B q − 1 \mathbf{B}_q^{-1} B q − 1 u q \mathbf{u}_q u q
Resolving the imposed field parallel and perpendicular to the pair axis gives the base field explicitly:
E n r = e z − cos Θ R ^ 1 + s E σ 21 δ + cos Θ R ^ 1 − 2 s E σ 21 δ . (179) \mathbf{E}_{\mathrm{nr}} =
\frac{\mathbf{e}_z-\cos\Theta\,\hat{\mathbf{R}}}{1+s_{\mathrm{E}}\sigma_{21}\delta}
+
\frac{\cos\Theta\,\hat{\mathbf{R}}}{1-2s_{\mathrm{E}}\sigma_{21}\delta}. \tag*{(179)} E nr = 1 + s E σ 21 δ e z − cos Θ R ^ + 1 − 2 s E σ 21 δ cos Θ R ^ . (179) F = sin 2 Θ ( 1 + s E σ 21 δ ) 2 + cos 2 Θ ( 1 − 2 s E σ 21 δ ) 2 , G = cos 2 Θ ( 1 − 2 s E σ 21 δ ) 2 . (180) F =
\frac{\sin^2\Theta}{(1+s_{\mathrm{E}}\sigma_{21}\delta)^2}
+
\frac{\cos^2\Theta}{(1-2s_{\mathrm{E}}\sigma_{21}\delta)^2},
\qquad
G = \frac{\cos^2\Theta}{(1-2s_{\mathrm{E}}\sigma_{21}\delta)^2}. \tag*{(180)} F = ( 1 + s E σ 21 δ ) 2 sin 2 Θ + ( 1 − 2 s E σ 21 δ ) 2 cos 2 Θ , G = ( 1 − 2 s E σ 21 δ ) 2 cos 2 Θ . (180) Expanding these formulas in δ \delta δ − 2 s E σ 21 ( 1 − 3 cos 2 Θ ) -2s_{\mathrm{E}}\sigma_{21}(1-3\cos^2\Theta) − 2 s E σ 21 ( 1 − 3 cos 2 Θ ) (98) and (99) of the main text. For s H = 0 s_{\mathrm{H}}=0 s H = 0 F F F
Crossover angle at finite separation
This calculation supplements Sec. VIII C and the kink in Figure 8 . It converts the local-field condition in Eq. (101) into the pair angle measured from the imposed field. For s H = 1 s_{\mathrm{H}}=1 s H = 1 X − X_- X − Y + Y_+ Y +
G F = 1 + δ / 2 3 ( 1 − δ / 2 ) . (181) \frac{G}{F}=\frac{1+\delta/2}{3(1-\delta/2)}. \tag*{(181)} F G = 3 ( 1 − δ /2 ) 1 + δ /2 . (181) Using Eq. (180) , the pair-axis angle at the crossing is
cos 2 Θ c r o s s = ( 1 + δ / 2 ) ( 1 − 2 s E σ 21 δ ) 2 ( 1 + δ / 2 ) ( 1 − 2 s E σ 21 δ ) 2 + 2 ( 1 − δ ) ( 1 + s E σ 21 δ ) 2 . (182) \cos^2\Theta_{\mathrm{cross}}=
\frac{(1+\delta/2)(1-2s_{\mathrm{E}}\sigma_{21}\delta)^2}
{(1+\delta/2)(1-2s_{\mathrm{E}}\sigma_{21}\delta)^2+2(1-\delta)(1+s_{\mathrm{E}}\sigma_{21}\delta)^2}. \tag*{(182)} cos 2 Θ cross = ( 1 + δ /2 ) ( 1 − 2 s E σ 21 δ ) 2 + 2 ( 1 − δ ) ( 1 + s E σ 21 δ ) 2 ( 1 + δ /2 ) ( 1 − 2 s E σ 21 δ ) 2 . (182) For small δ \delta δ
cos 2 Θ c r o s s = 1 3 + δ 3 − 4 3 s E σ 21 δ + O ( δ 2 ) . (183) \cos^2\Theta_{\mathrm{cross}}=\frac{1}{3}+\frac{\delta}{3}-\frac{4}{3}s_{\mathrm{E}}\sigma_{21}\delta+\mathcal{O}(\delta^2). \tag*{(183)} cos 2 Θ cross = 3 1 + 3 δ − 3 4 s E σ 21 δ + O ( δ 2 ) . (183) The term δ / 3 \delta/3 δ /3 − 4 s E σ 21 δ / 3 -4s_{\mathrm{E}}\sigma_{21}\delta/3 − 4 s E σ 21 δ /3 σ 21 = − 1 / 2 \sigma_{21}=-1/2 σ 21 = − 1/2
Θ c r o s s = Θ m − 3 δ 2 2 + O ( δ 2 ) , Θ c r o s s H = Θ m − δ 2 2 + O ( δ 2 ) . (184) \Theta_{\mathrm{cross}}=\Theta_{\mathrm{m}}-\frac{3\delta}{2\sqrt{2}}+\mathcal{O}(\delta^2),\qquad
\Theta_{\mathrm{cross}}^{\mathrm{H}}=\Theta_{\mathrm{m}}-\frac{\delta}{2\sqrt{2}}+\mathcal{O}(\delta^2). \tag*{(184)} Θ cross = Θ m − 2 2 3 δ + O ( δ 2 ) , Θ cross H = Θ m − 2 2 δ + O ( δ 2 ) . (184) These finite-separation displacements of the unexpanded crossover approach zero with R − 3 R^{-3} R − 3 O ( R − 6 ) \mathcal{O}(R^{-6}) O ( R − 6 )
Regular asymptotic steady angular velocities away from onset
This subsection derives the angular-speed corrections in Eqs. (104) and (105) of Sec. IX, and their overlap with the near-onset theory. This calculation uses the regular asymptotic expansion away from onset, at fixed λ 0 > 0 \lambda_0>0 λ 0 > 0 δ / λ 0 ≪ 1 \delta/\lambda_0 \ll1 δ / λ 0 ≪ 1 P ⊥ ( 0 ) = M a λ 0 P_{\perp}^{(0)}=\mathrm{Ma}\sqrt{\lambda_0} P ⊥ ( 0 ) = Ma λ 0 P i ( 0 ) = P ⊥ ( 0 ) e ρ i + ( ϵ 21 − M a ) e z P_i^{(0)}=P_{\perp}^{(0)}e_{\rho i}+(\epsilon_{21}-\mathrm{Ma})e_z P i ( 0 ) = P ⊥ ( 0 ) e ρ i + ( ϵ 21 − Ma ) e z G i , ϕ = 0 G_{i,\phi}=0 G i , ϕ = 0
v i = Π ⋅ P j ( 0 ) , w i = Π ⋅ e ϕ j , (185) v_i=\boldsymbol{\Pi}\cdot P_j^{(0)},\qquad w_i=\boldsymbol{\Pi}\cdot e_{\phi j}, \tag*{(185)} v i = Π ⋅ P j ( 0 ) , w i = Π ⋅ e ϕ j , (185) and use ρ , ϕ , z \rho,\phi,z ρ , ϕ , z M a − 1 [ − P i ( 0 ) × v i + ( P ⊥ ( 0 ) / 2 ) w i ] \mathrm{Ma}^{-1}[-P_i^{(0)}\times v_i+(P_{\perp}^{(0)}/2)w_i] Ma − 1 [ − P i ( 0 ) × v i + ( P ⊥ ( 0 ) /2 ) w i ]
P i , ρ ( 1 ) = ( M a − ϵ 21 ) v i , ρ − M a 2 P ⊥ ( 0 ) v i , z + ( P ⊥ ( 0 ) ) 2 − M a 2 4 P ⊥ ( 0 ) w i , ϕ , (186) P_{i,\rho}^{(1)}=(\mathrm{Ma}-\epsilon_{21})v_{i,\rho}-\frac{\mathrm{Ma}^{2}}{P_{\perp}^{(0)}}v_{i,z}+\frac{(P_{\perp}^{(0)})^{2}-\mathrm{Ma}^{2}}{4P_{\perp}^{(0)}}w_{i,\phi}, \tag*{(186)} P i , ρ ( 1 ) = ( Ma − ϵ 21 ) v i , ρ − P ⊥ ( 0 ) Ma 2 v i , z + 4 P ⊥ ( 0 ) ( P ⊥ ( 0 ) ) 2 − Ma 2 w i , ϕ , (186) P i , z ( 1 ) = P ⊥ ( 0 ) v i , ρ − ( M a + ϵ 21 ) v i , z − M a 2 w i , ϕ . (187) P_{i,z}^{(1)}=P_{\perp}^{(0)}v_{i,\rho}-(\mathrm{Ma}+\epsilon_{21})v_{i,z}-\frac{\mathrm{Ma}}{2}w_{i,\phi}. \tag*{(187)} P i , z ( 1 ) = P ⊥ ( 0 ) v i , ρ − ( Ma + ϵ 21 ) v i , z − 2 Ma w i , ϕ . (187) The phase coordinate is chosen so that P i , ϕ ( 1 ) = 0 P_{i,\phi}^{(1)}=0 P i , ϕ ( 1 ) = 0 X − X_- X − Y + Y_+ Y + Θ = π / 2 \Theta=\pi/2 Θ = π /2 ∣ Ω i ∣ = λ 0 + δ ω i ( 1 ) + O ( δ 2 ) |\Omega_i|=\sqrt{\lambda_0}+\delta\omega_i^{(1)}+O(\delta^2) ∣ Ω i ∣ = λ 0 + δ ω i ( 1 ) + O ( δ 2 )
ω i ( 1 ) = − ( P ⊥ ( 0 ) ) 2 + M a 2 P ⊥ ( 0 ) M a ( v i , z + 1 4 w i , ϕ ) . (188) \omega_i^{(1)}=-\frac{(P_{\perp}^{(0)})^2+\mathrm{Ma}^2}{P_{\perp}^{(0)}\mathrm{Ma}}\left(v_{i,z}+\frac{1}{4}w_{i,\phi}\right). \tag*{(188)} ω i ( 1 ) = − P ⊥ ( 0 ) Ma ( P ⊥ ( 0 ) ) 2 + Ma 2 ( v i , z + 4 1 w i , ϕ ) . (188) This identity follows by combining the radial dipole contribution P i , ρ ( 1 ) / M a P_{i,\rho}^{(1)}/\mathrm{Ma} P i , ρ ( 1 ) / Ma v i , ρ v_{i,\rho} v i , ρ
For X − X_- X − e ρ 1 = e x e_{\rho1}=e_x e ρ 1 = e x e ρ 2 = − e x e_{\rho2}=-e_x e ρ 2 = − e x w i , ϕ = − 1 w_{i,\phi}=-1 w i , ϕ = − 1
1 2 ( v 1 , z + v 2 , z ) = ( ϵ 21 − M a ) ( 1 − 3 cos 2 Θ ) . (189) \frac{1}{2}(v_{1,z}+v_{2,z})=(\epsilon_{21}-\mathrm{Ma})(1-3\cos^2\Theta). \tag*{(189)} 2 1 ( v 1 , z + v 2 , z ) = ( ϵ 21 − Ma ) ( 1 − 3 cos 2 Θ ) . (189) Averaging the two positive speeds gives Ω ‾ X − \overline{\Omega}_{X_-} Ω X − (104) of the main text. For Y + Y_+ Y + Θ = π / 2 \Theta=\pi/2 Θ = π /2 e ρ i = e y e_{\rho i}=e_y e ρ i = e y w i , ϕ = − 2 w_{i,\phi}=-2 w i , ϕ = − 2 v i , z = ϵ 21 − M a v_{i,z}=\epsilon_{21}-\mathrm{Ma} v i , z = ϵ 21 − Ma Y + Y_+ Y + (105) of the main text. Oblique X − X_- X − sin Θ cos Θ \sin\Theta\cos\Theta sin Θ cos Θ
The overlap with the near-onset theory is obtained by taking δ ≪ λ 0 ≪ 1 \delta\ll\lambda_0\ll1 δ ≪ λ 0 ≪ 1 M a = M a c / ( 1 + λ 0 ) \mathrm{Ma}=\mathrm{Ma}_c/(1+\lambda_0) Ma = Ma c / ( 1 + λ 0 )
Ω ‾ k = λ 0 + δ m k 2 λ 0 + O ( δ λ 0 ) + O ( δ 2 λ 0 − 3 / 2 ) , k = X − , Y + . (190) \overline{\Omega}_k=\sqrt{\lambda_0}+\frac{\delta m_k}{2\sqrt{\lambda_0}}+O(\delta\sqrt{\lambda_0})+O(\delta^2\lambda_0^{-3/2}),\qquad k=X_-,Y_+. \tag*{(190)} Ω k = λ 0 + 2 λ 0 δ m k + O ( δ λ 0 ) + O ( δ 2 λ 0 − 3/2 ) , k = X − , Y + . (190) The same expression follows by expanding λ 0 + δ m k \sqrt{\lambda_0+\delta m_k} λ 0 + δ m k λ 0 \lambda_0 λ 0 λ 0 = O ( δ ) \lambda_0=O(\delta) λ 0 = O ( δ )
