2. The swimmer and equations of motion

We consider a recently proposed swimmer, the Quadroar (Jalali, Alam & Mousavi Reference Jalali, Alam and Mousavi2014), which is composed of four rotating disks of radius $a$ at the ends of its two axles, and a chassis of variable length. Each axle has a length of $2b$ and the length of the chassis $2l+2s(t)$ is varied by a linear actuator (figure 1 a). The chassis and axles make an I-shape planar frame. A body-fixed coordinate system $(x_{1},x_{2},x_{3})$ with the unit base vectors $(\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3})$ is chosen to describe the kinematics and dynamics of the swimmer’s translational and rotational movements. The $x_{1}$ -axis is along the chassis, the $x_{2}$ -axis is parallel to the axles and the $x_{3}$ -axis is normal to the plane of the I-frame. The relative orientation of the $n$ th disk with respect to the swimmer’s frame is measured by the angle $-{\rm\pi}\leqslant {\it\vartheta}_{n}\leqslant +{\rm\pi}$ that the plane of the disk makes with the $(x_{1},x_{2})$ coordinate plane. We define the global coordinate frame $(X_{1},X_{2},X_{3})$ with the unit base vectors $(\boldsymbol{E}_{1},\boldsymbol{E}_{2},\boldsymbol{E}_{3})$ , and denote the position vector of the swimmer’s centre of mass by $\boldsymbol{X}_{c}(t)$ . The translational velocity of the swimmer thus becomes $\boldsymbol{v}_{c}=\dot{\boldsymbol{X}_{c}}$ . Here, an overdot stands for $\text{d}/\text{d}t$ . We will represent all physical quantities in terms of the unit vectors $\boldsymbol{e}_{i}$ in the body-fixed coordinate frame, and will explicitly state if we switch to the global coordinate system.

Figure 1. (a) The geometry of the Quadroar. The centre of the body-fixed coordinate system $(x_{1},x_{2},x_{3})$ is at the centre of mass of the swimmer. (b) The geometry of the swimmer as seen along the positive $x_{2}$ -axis, which is parallel to the axles. The chassis performs its open $\rightleftarrows$ close function according to a harmonic function with a constant frequency ${\it\omega}_{s}$ . The front and rear disks counter-rotate with the base angular velocities $\dot{{\it\vartheta}}_{1}=\dot{{\it\vartheta}}_{2}={\it\omega}_{s}/2$ and $\dot{{\it\vartheta}}_{3}=\dot{{\it\vartheta}}_{4}=-{\it\omega}_{s}/2$ , respectively, but their rotational frequencies can be shifted by ${\rm\Delta}{\it\nu}_{n}$ .

The main difference between the Quadroar and other swimmer designs such as those that use linked spheres is the point torque that each rotating disk of the Quadroar exerts on the background fluid. Consequently, the streaming of the background fluid in response to the combined translational and rotational motions of the disks (propellers) is a combination of Stokeslets and rotlets. Jalali et al. (Reference Jalali, Alam and Mousavi2014) ignored the hydrodynamic interactions between the disks of the swimmer by assuming that they are sufficiently far from each other. Since this study is going to address the near and far flow fields generated by the swimmer, hydrodynamic interactions of the four disks cannot be ignored. We thus derive new equations of motion for the swimming of the Quadroar taking full account of hydrodynamic interactions. Nonetheless, we retain one of our basic simplifying assumptions: each disk interacts with the background fluid through a point force and a point torque, both applied at the centre of the disk. Therefore, in our analysis, the disks still need to be far from each other such that geometrical effects on streamlines due to the ‘finite sizes’ of the disks can be ignored.

The Quadroar can propel both in stepwise and continuous operation modes. In this study, we consider continuous operation modes that result in two-dimensional orbits of the swimmer in its sagittal $(x_{1},x_{3})$ -plane. In the operation mode that we are interested in, the chassis expands and contracts according to the harmonic function $s(t)=s_{0}[1-\cos ({\it\omega}_{s}t)]/2$ , and the front and rear disks counter-rotate with the base frequency ${\it\omega}_{s}/2$ (figure 1 b). Our control input is the small frequency shift ${\rm\Delta}{\it\nu}$ that detunes the rotational speeds of the front and rear disks as

The swimming dynamics is therefore associated with two time scales: the fast time $T_{\mathit{fast}}=2{\rm\pi}/{\it\omega}_{s}$ and the slow time $T_{\mathit{slow}}=2{\rm\pi}/{\rm\Delta}{\it\nu}$ . Jalali et al. (Reference Jalali, Alam and Mousavi2014) showed that the Quadroar strokes rectilinearly along its $x_{3}$ body axis for ${\rm\Delta}{\it\nu}=0$ , and swims on quasi-periodic orbits for non-zero frequency shifts. While depending on the geometrical aspects of the Quadroar details of trajectories maybe different, quasi-periodic orbits also exist when hydrodynamic interactions between disks are included. Quasi-periodic orbits have remarkable consequences on the streaming of the background fluid, which we explore here.

The relative position vector of the $n$ th disk with respect to the centre of mass of the swimmer is defined by $\boldsymbol{r}_{n}$ . We have

Due to the expanding and contracting motion of the chassis (body link), each disk of the Quadroar has a relative velocity $\boldsymbol{v}_{\mathit{rel},n}$ with respect to the centre of mass. The relative velocities are computed by taking the time derivatives of $\boldsymbol{r}_{n}$ in the body-fixed coordinate system. For the front disks $(n=1,2)$ and rear disks $(n=3,4)$ one has $\boldsymbol{v}_{\mathit{rel},1}=\boldsymbol{v}_{\mathit{rel},2}={\dot{s}}\boldsymbol{e}_{1}$ and $\boldsymbol{v}_{\mathit{rel},3}=\boldsymbol{v}_{\mathit{rel},4}=-{\dot{s}}\boldsymbol{e}_{1}$ , respectively. Here ${\dot{s}}$ is the expansion or contraction rate of the chassis.

The angular velocity of the swimmer, which is also the angular velocity of the body-fixed coordinate system, is denoted by ${\bf\omega}_{\mathit{body}}(t)$ . It is related to the 1-2-3 sequence of Euler angles ${\bf\alpha}=({\it\phi},{\it\theta},{\it\psi})$ through the following transformation

where $s_{{\it\gamma}}$ and $c_{{\it\gamma}}$ denote $\sin ({\it\gamma})$ and $\cos ({\it\gamma})$ , respectively. The transformation between body-fixed and global coordinate systems is carried out using the rotation matrix

For instance, the velocity of the centre of mass in the body-fixed coordinate system reads $\unicode[STIX]{x1D64D}\boldsymbol{\cdot }\dot{\boldsymbol{X}_{c}}$ . The absolute linear and angular velocities of the $n$ th disk are thus determined from

We define the isotropic tensor $\unicode[STIX]{x1D642}=(32/3)a^{3}\unicode[STIX]{x1D644}$ , with $\unicode[STIX]{x1D644}$ being the identity matrix and the translation matrix (Jalali et al. Reference Jalali, Alam and Mousavi2014)

At its hydrodynamic centre, the $n$ th disk exerts the force vector $\boldsymbol{f}_{n}$ and torque ${\bf\tau}_{n}$ on the fluid. The reactions of these are the drag forces and torques. Within a non-inertial streaming fluid, the following analytic expressions are known (Happel & Brenner Reference Happel and Brenner1983)

where $\boldsymbol{u}(\boldsymbol{X},t)$ and $2{\it\bf\Omega}(\boldsymbol{X},t)=\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}$ are the streaming and vorticity fields of the background fluid with the dynamic viscosity ${\it\mu}$ , and we have

It is noted that the drag forces and torques occur when the relative linear and angular velocities within the parentheses in (2.8a,b ) are non-zero. The effects of $\boldsymbol{u}_{n}$ and ${\it\bf\Omega}_{n}$ cannot be neglected for compact swimmers whose own propellers are close to each other and have hydrodynamic interactions.

For a neutrally buoyant swimmer in low Reynolds number conditions, $\mathit{Re}\rightarrow 0$ , the drag force and torque vanish and we obtain

These equations are still incomplete for the determination of the translational and rotational dynamics of the swimmer as they contain the unknown velocities $\boldsymbol{u}_{n}$ and spins ${\it\bf\Omega}_{n}$ that are developed in the background fluid because of the activity of the swimmer itself: each disk generates a flow and influences the operation of other propellers. If we assume that the $k$ th disk exerts a point force and torque (at its hydrodynamic centre) on the fluid, the velocity field that it generates at $\boldsymbol{X}$ will be the combination of a Stokeslet and a rotlet as (e.g. Lopez & Lauga Reference Lopez and Lauga2014)

Operating (2.12) with $\boldsymbol{{\rm\nabla}}\times$ gives the vorticity field:

The velocity and vorticity fields induced at the hydrodynamic centre of the $n$ th disk by three other disks are thus determined from

where $\boldsymbol{X}_{kn}=\boldsymbol{r}_{n}-\boldsymbol{r}_{k}$ and $z_{kn}=(\boldsymbol{X}_{kn}\boldsymbol{\cdot }\boldsymbol{X}_{kn})^{1/2}$ . In these equations, the forces $\boldsymbol{f}_{k}$ and torques ${\bf\tau}_{k}$ linearly depend on $\boldsymbol{u}_{k}$ and ${\it\bf\Omega}_{k}$ through relations in (2.8a,b ). We collect the 30 unknown components of the vectors $\boldsymbol{v}_{c}$ , ${\bf\omega}_{\mathit{body}}$ , $\boldsymbol{u}_{n}$ and ${\it\bf\Omega}_{n}\;(n=1,2,3,4)$ in the vector $\boldsymbol{w}$ , and combine (2.14) and (2.15) with (2.10) and (2.11) to obtain a system of linear equations

in terms of $\boldsymbol{w}$ and $t$ , with $\boldsymbol{L}$ being a vectorial function of dimension $30\times 1$ . We then analytically calculate the Jacobian of $\boldsymbol{L}$ and transform (2.16) to

Here the matrix $\unicode[STIX]{x1D63C}(t)$ and the right-hand side vector $\boldsymbol{d}(t)$ explicitly depend on $t$ through the control variables $s(t)$ and ${\it\vartheta}_{n}(t)$ . At each time step, the linear system of (2.17) is solved using a standard LU decomposition method with pivoting. After calculating $\boldsymbol{v}_{c}$ and ${\bf\omega}_{\mathit{body}}$ , the position vector $\boldsymbol{X}_{c}=X_{c,i}\boldsymbol{E}_{i}$ and orientation angles of the swimmer are found by the direct numerical integration of

with the superscript T denoting transpose. We now express all physical quantities in terms of the unit vectors $\boldsymbol{E}_{i}$ of the global coordinate system. This is simply done by left-multiplying all vectors in the $(x_{1},x_{2},x_{3})$ frame by $\unicode[STIX]{x1D64D}^{\text{T}}$ . The flow field generated by the Quadroar is a superposition of the Stokeslets and rotlets of the four disks:

where $\boldsymbol{f}_{n}$ and ${\bf\tau}_{n}$ are computed using the components of $\boldsymbol{w}(t)$ . We study the flow characteristics by following the Lagrangian trajectories of passive tracers. The equation of motion for a tracer becomes

which is integrated in the global coordinate system simultaneous with the swimmer’s equations of motion in (2.18a,b ). Due to the coupled rotational and translational motions of the swimmer and their effect on the background fluid and tracers, the equations of motion are highly nonlinear and we need an accurate integrator to avoid spurious mixing due to numerical errors. Leap frog integrators devised for second-order ordinary differential equations do not work well here because we do not have evolution equations for the velocities – there is no acceleration in the low Reynolds flow regimes. We therefore adopt the RK78 integrator (Fehlberg Reference Fehlberg1968) that can keep relative integration errors at $O(10^{-14})$ over long time scales that a quasi-periodic orbit becomes dense in its invariant manifold or chaotic orbits uniformly sample their invariant measure.

3. Flow field

The Quadroar is a three-dimensional swimmer with rich orbital structure in the phase space. Nonetheless, to study the quantitative and qualitative effects of swimming on environmental fluid, and gain insight into the regions of influence of microswimmers, we investigate two-dimensional quasi-periodic motions in the body-fixed $(x_{1},x_{3})$ coordinate system. The parameters of our model swimmer have been set to $a=1$ , $b=4$ , $l=4$ , ${\it\omega}_{s}=1$ , and $s_{0}/l=1/2$ . We assume that all disks start their rotations in phase, ${\it\vartheta}_{n}=0\;(n=1,2,3,4)$ , and the initial conditions for the position and orientation of the swimmer are set to $\boldsymbol{X}_{c}=\mathbf{0}$ and ${\bf\alpha}=\mathbf{0}$ . Depending on our choice of the control parameter ${\rm\Delta}{\it\nu}$ in (2.1a,b ), the swimmer can take different courses of quasi-periodic orbits.

The simplest motion is rectilinear with ${\rm\Delta}{\it\nu}=0$ : the front and rear disks counter-rotate as the chassis periodically expands and contracts. Disks simultaneously become parallel and perpendicular to the swimmer’s body plane $(x_{1},x_{2})$ at the times $t=j\;T_{\mathit{fast}}$ and $t=(2j+1)T_{\mathit{fast}}/2$ ( $j=0,1,2,\ldots$ ), respectively. When the chassis is expanding, ${\dot{s}}>0$ , the inclination angles of the front and rear disks (with respect to the swimmer’s body plane) vary in the intervals $0\leqslant {\it\vartheta}_{1}={\it\vartheta}_{2}\leqslant {\rm\pi}/2$ and $-{\rm\pi}/2\leqslant {\it\vartheta}_{3}={\it\vartheta}_{4}\leqslant 0$ , respectively. In such conditions, the swimmer has the conformation of figure 2(a). During the contraction phase of the chassis, ${\dot{s}}<0$ , the front and rear disks respectively evolve in the intervals ${\rm\pi}/2\leqslant {\it\vartheta}_{1}={\it\vartheta}_{2}\leqslant {\rm\pi}$ and $-{\rm\pi}\leqslant {\it\vartheta}_{3}={\it\vartheta}_{4}\leqslant -{\rm\pi}/2$ and the swimmer has the conformation of figure 2(b). We have used (2.19) to compute the flow field $\boldsymbol{u}(\boldsymbol{X},t)$ induced by the swimmer. Streamlines and velocity vectors have been displayed in figure 3 for several values of $0<t<T_{\mathit{fast}}$ during a full cycle of the swimmer’s stroke. Since there is no net motion/flow in the $x_{2}$ direction, we have visualised $\boldsymbol{u}$ only in the $(X_{1},X_{3})$ slice of the configuration space that overlaps with the $(x_{1},x_{3})$ coordinate plane. Another justification for selecting the $(X_{1},X_{3})$ plane is to avoid the geometrical effects of the disks, which have finite sizes: our methodology of modelling the disks by $\boldsymbol{f}_{k}$ and ${\bf\tau}_{k}$ is supposed to give reasonable results for $\boldsymbol{u}$ if $|\boldsymbol{X}-\boldsymbol{X}_{c}-\boldsymbol{r}_{n}|\gg a$ . This condition is fulfilled in the $(x_{1},x_{3})$ plane.

Figure 2. Side views of the Quadroar swimmer during a rectilinear motion: the front and rear disks counter-rotate with the same angular velocities, and their inclinations with respect to the chassis are identical but with different signs. (a) and (b) respectively show the swimmer’s conformations during the expansion ( ${\dot{s}}>0$ ) and contraction ( ${\dot{s}}<0$ ) phases of the chassis.

Figure 3. Snapshots of the streamlines of the flow field $\boldsymbol{u}(\boldsymbol{X},t)$ induced by the Quadroar swimmer that starts its motion form $(X_{1},X_{3})=(0,0)$ and strokes along the negative $X_{3}$ -axis. The chassis of the swimmer with variable length $2l+2s(t)$ has been shown by a thick red bar. For each panel, the conformation of the swimmer’s disks can be determined in terms of $\text{sign}({\dot{s}})$ and figure 2. The inclinations of the disks with respect to the swimmer’s chassis can be computed using the following equations: ${\it\vartheta}_{1}={\it\vartheta}_{2}={\it\omega}_{s}t/2$ and ${\it\vartheta}_{3}={\it\vartheta}_{4}=-{\it\omega}_{s}t/2$ . (a)  $t=(1/64)T_{\mathit{fast}}$ , ${\dot{s}}>0$ , (b)  $t=(8/64)T_{\mathit{fast}}$ , ${\dot{s}}>0$ , (c)  $t=(16/64)T_{\mathit{fast}}$ , ${\dot{s}}>0$ , (d)  $t=(32/64)T_{\mathit{fast}}$ , ${\dot{s}}=0$ , (e)  $t=(48/64)T_{\mathit{fast}}$ , ${\dot{s}}<0$ , (f)  $t=(63/64)T_{\mathit{fast}}$ , ${\dot{s}}<0$ .

Each stroke cycle begins with the expansion of the chassis and two lateral vortices are formed on the anterior side of the swimmer. These two vortices are pushed forward and shrink in size as the inclination angles ${\it\vartheta}_{n}$ increase (figure 3 a,b). For $|{\it\vartheta}_{n}|\approx {\rm\pi}/4$ , the lateral vortices disappear and the flow field takes a hyperbolic structure where the fluid is pumped in towards the centre of mass of the swimmer along the $X_{3}$ -axis, and is pumped out along the $X_{1}$ -axis (figure 3 c). During this transitional phase, the swimmer behaves like a puller dipole swimmer (e.g. Elgeti, Winkler & Gompper Reference Elgeti, Winkler and Gompper2015). As the disks make right angles with respect to the chassis, two counter-rotating side vortices are generated by the disks (figure 3 d). When the chassis contracts, the swimmer behaves like a pusher dipole: the fluid is pumped in along the $X_{1}$ -axis and pumped out along the $X_{3}$ -axis as in pushers (figure 3 e). The flow structure is hyperbolic for $|{\it\vartheta}_{n}|\approx {\rm\pi}/4$ , and lateral vortices occur behind the swimmer for smaller inclinations of the disks (figure 3 f). The Quadroar is therefore a hybrid swimmer that inherits the features of both pusher and puller dipole swimmers. The streamlines displayed in figure 3 resemble the oscillatory flow field of a single-celled swimming alga called C. reinhardtii (Guasto, Johnson & Gollub Reference Guasto, Johnson and Gollub2010).

Figure 4. (a) The quasi-periodic orbit of the swimmer integrated from $(X_{1},X_{3})=(0,0)$ over the time interval $0\leqslant t\leqslant t_{max}=6.25T_{\mathit{slow}}$ . The orbit precesses and becomes dense in an annular region. (b) A close-up snapshot of the orbit near the swimmer’s position at $t_{max}$ . The orientation of the body-fixed coordinate frame $(x_{1},x_{3})$ is displayed at $t_{max}$ . (c) The streamlines and local velocity vectors of the flow field $\boldsymbol{u}(\boldsymbol{X},t)$ at $t=6.25T_{\mathit{slow}}-(7/8)T_{\mathit{fast}}$ . (d) Same as (c) but for $t=6.25T_{\mathit{slow}}-(3/4)T_{\mathit{fast}}$ . The thick bar shows the chassis of the swimmer.

4. Chaotic mixing

The rectilinear mode of swimming helps us gain insight into the flow physics around microswimmers, but it is not useful for mixing due to two obvious reasons: (i) a rectilinear trajectory has measure zero in the two-dimensional configuration space, and (ii) fluid particles affected by a distant swimmer passing on a straight line move on loop-like paths (Pushkin & Yeomans Reference Pushkin and Yeomans2013), which are not effective in the stirring of environmental fluid. We therefore turn our attention to curved quasi-periodic orbits of non-zero measure. We set ${\rm\Delta}{\it\nu}=0.00125$ that generates the precessing orbit of figure 4(a) for the swimmer’s centre of mass. The trajectory has inverted multi-loop turns aligned towards the centre of the orbit and becomes dense in an annular region $\mathscr{A}$ , which is the projection of a three-dimensional torus $\mathscr{M}$ in the $(X_{1},X_{3},{\it\theta})$ -space on the $(X_{1},X_{3})$ -plane. In the literature of dynamical systems theory, $\mathscr{M}$ is called the ‘invariant torus’ of a quasi-periodic orbit. We have integrated the equations of motion from $t=0$ to $t_{max}=6.25T_{\mathit{slow}}$ . Figure 4(b) shows a close-up of the trajectory and highlights it over the time interval $6.03T_{\mathit{slow}}\leqslant t\leqslant t_{max}$ . We have also shown in this figure the location and orientation of the swimmer’s chassis and its body-fixed coordinate system at $t_{max}$ .

Near the state of the swimmer shown in figure 4(b), we have plotted the streamlines and local velocity vectors at the two instants $t=6.25T_{\mathit{slow}}-(7/8)T_{\mathit{fast}}$ (figure 4 c) and $t=6.25T_{\mathit{slow}}-(3/4)T_{\mathit{fast}}$ (figure 4 d). We have only shown the velocity field in the $(x_{1},x_{3})$ body plane that passes through the centre of mass of the swimmer and overlaps with the global $(X_{1},X_{3})$ plane. It can be seen how the swimmer pumps the fluid while making a spiral turn. Here the magnitudes of ${\it\vartheta}_{1,2}$ and ${\it\vartheta}_{3,4}$ differ due to gradual phase shift generated by ${\rm\Delta}{\it\nu}$ , and streamlines are rapidly evolving. The twisted hyperbolic structure of figure 4(c) occurs when the swimmer strokes like a puller dipole, and four separatrices intersect at a saddle node which is close to the swimmer’s centre of mass. The flow stream of figure 4(d) originates around one axle and sinks towards the other one. This pattern has no analog in rectilinear motion and is a consequence of phase shift between front and rear disks and spiralling motion. These patterns are few demonstrative examples which show how a swimmer moving on a curved path induces complex, evolving topologies of streamlines. Depending on the phase shift ${\it\vartheta}_{1,2}+{\it\vartheta}_{3,4}={\rm\Delta}{\it\nu}\;t$ , the swimmer’s trajectory and its orientation, new topologies can emerge.

Figure 5. Dispersion of a packet of passive tracers by a single quasi-periodic orbit swimmer. The packet has been split to four subdomains of different colours. (a) The initial condition of the packet, and its subsequent states at $t/T_{\mathit{slow}}=3.75$ , 6.25 and 10 where $T_{\mathit{slow}}=2{\rm\pi}/{\rm\Delta}{\it\nu}$ . The packet is folded and stretched in the flow field of the swimmer of figure 4. (b) A highly folded and stretched intermediary state of the patch at $t=22.5T_{\mathit{slow}}$ . Contact lines between four subdomains are being destroyed: (c)  $t=50T_{\mathit{slow}}$ . The contact lines between the four subdomains of the initial patch have been completely destroyed and the mixing process has started. (d) The fully mixed state of the patch at $t=125T_{\mathit{slow}}$ . Tracers have been dispersed uniformly within the annular region $\mathscr{A}$ filled by the swimmer’s orbit. The initial packet is also superimposed in panels (c) and (d) for comparison. Note the different scales of panels.

To understand the mixing process and its efficiency, we follow the Lagrangian trajectories of a packet of tracers as the orbit of the swimmer illustrated in figure 4(a) becomes dense in the annular region $\mathscr{A}$ . We choose a uniformly distributed set of 2500 passive tracers in the rectangular region $\mathscr{D}_{0}=\{(X_{1},X_{3})|-70\leqslant X_{1}\leqslant -67.5,-1.25\leqslant X_{3}\leqslant 1.25\}$ , which lies well in the ring of the swimmer’s orbit. To visualise the mixing process and deformation of fluid elements, we split $\mathscr{D}_{0}$ to four equal square subdomains and tag their corresponding tracers in a clockwise order by filled blue, green, black and red dots (figure 5 a). We then release the swimmer from $(X_{1},X_{3})=(0,0)$ and integrate (2.18a,b ) for the swimmer and (2.21) for the tracers. The fluid element $\mathscr{D}_{0}$ starts to deform and is mapped to the set $\mathscr{D}(t)$ as the time elapses. Figure 5(a) demonstrates the shape of $\mathscr{D}(t)$ , with $\mathscr{D}(0)=\mathscr{D}_{0}$ , at four snapshots in time. It is seen that the fluid element and the contact lines between its four subdomains are deformed gradually from $t=0$ to $t=6.25T_{\mathit{slow}}$ but the element keeps its integrity.

The element is highly stretched and folded later at $t=10T_{\mathit{slow}}$ and starts to develop tails. Over this period, the swimmer has completed approximately 21 turning loops of its quasi-periodic orbit. The folded patch evolves to thin wavy structures as seen in figure 5(b) and the contact lines between the four subdomains are destroyed, allowing tracers to gradually distance from their associates. The set $\mathscr{D}(t)$ shown in figure 5(c) at $t=50T_{\mathit{slow}}$ no longer exhibits the properties of a manifold and is topologically similar to strange attractors of dissipative dynamical systems, although it is just an intermediary and not a final stage of an ongoing mixing process. The swimmer completely disperses tracers in $\mathscr{A}$ as $t\rightarrow \infty$ . The tracers of the four subdomains of $\mathscr{D}_{0}$ have been well mixed by $t=125T_{\mathit{slow}}$ (figure 5 d).

Figure 6. (a) The profile of $T_{\mathit{slow}}{\it\sigma}_{max}(\mathscr{L}_{0},50T_{\mathit{slow}})$ along the line element $\mathscr{L}_{0}$ . (b) The state of the line element represented by 1000 tracer particles at $t=0$ (blue solid line) and $t=50T_{\mathit{slow}}$ (black dots). The element has been disrupted along the path of the swimmer (cf. figure 4 a).

Various measures can be utilized to quantify the mixing process, including entropy, Lyapunov exponents and the mix-norm of advected scalar fields by chaotic flows (Mathew, Mezić & Petzold Reference Mathew, Mezić and Petzold2005). In this study, we choose the line element $\mathscr{L}_{0}=\{(X_{1},X_{3})|X_{3}=0,-100\leqslant X_{1}\leqslant -40\}$ , and compute its corresponding finite-time largest Lyapunov exponent, ${\it\sigma}_{max}$ , following the procedure of Chabreyrie, Chandre & Aubry (Reference Chabreyrie, Chandre and Aubry2011). A positive Lyapunov exponent implies chaos and exponential divergence of neighbouring tracers that lie on $\mathscr{L}_{0}$ . To compute Lyapunov exponents corresponding to a particle initially located at $\boldsymbol{X}_{0}=\boldsymbol{X}(0)$ , we integrate

along with (2.21) and compute the $3\times 3$ state transition matrix $\unicode[STIX]{x1D731}(t)$ and its three eigenvalues ${\it\lambda}_{k}\;(k=1,2,3)$ at $t=T=50T_{\mathit{slow}}$ . The largest Lyapunov exponent is calculated from ${\it\sigma}_{max}=T^{-1}max({\it\lambda}_{1},{\it\lambda}_{2},{\it\lambda}_{3})$ . Figure 6(a) shows the variation of $T_{\mathit{slow}}{\it\sigma}_{max}(\mathscr{L}_{0},T)$ . Results do not considerably change by increasing $T$ . We have also shown the deformed and disrupted state of $\mathscr{L}_{0}$ at $t=T$ . Our numerical experiments show that ${\it\sigma}_{max}$ is positive wherever the line element is stretched, and becomes minimum at the tips of the deformed lobes where tracer particles move together and the line element experiences the least stretching and maximum folding. The magnitude of the Lyapunov exponent is small because the swimmer spends a large fraction of time far from a test tracer: a quasi-periodic orbit is characterised by a short period needed to complete a full multi-loop tun and a subsequent long-range stroke, and a long period associated with its orbital precession. During each swimmer–tracer encounter, tracer particles are displaced but they remain almost stationary until the swimmer revisits that part of the configuration space. Since the precession rate of quasi-periodic orbits is ${\it\varpi}\sim O(T_{\mathit{slow}}^{-1})$ , the recurrence time for swimmer–tracer interaction is ${\sim}2{\rm\pi}/{\it\varpi}$ , which is long. Therefore, Lyapunov exponents are scaled by $T_{\mathit{slow}}^{-1}$ , and remain almost zero over the hibernation (stationary) phase of tracer particles. These two effects decrease the magnitude of ${\it\sigma}_{max}$ and increase the time scale of mixing substantially. The mixing speed can be enhanced using more swimmers on a given orbit, but hydrodynamic interactions between swimmers must be taken into account.

We have repeated our simulations for different choices of the frequency shift ${\rm\Delta}{\it\nu}$ . Our numerical experiments show that chaotic mixing occurs within the invariant torus of all quasi-periodic orbits, and the thickness of the chaotic layer correlates with the structure of the turning loops: the higher the number of loops in a turning multi-loop bundle, the greater the number of foldings that a test patch experiences, and consequently the wider the chaotic layer. The number of loops nonlinearly depends on the geometric parameters of the swimmer and ${\rm\Delta}{\it\nu}$ . The turning loop of the simplest quasi-periodic orbit in our library is a cardioid corresponding to ${\rm\Delta}{\it\nu}=0.00625$ (figure 7 a). For this swimming mode, the chaotic set forms a rotating chain-like structure in the configuration space. We have not observed any sign of chaos for periodic orbits. For ${\rm\Delta}{\it\nu}=0.005$ , the swimmer’s orbit is periodic with inverted single cardioids as turning loops (figure 7 b). It is seen that fluid elements are stretched and carried by the swimmer on a curved path, but they are not subject to successive foldings and do not diffuse in the annular region covered by the swimmer’s orbit. Full exploration of chaotic zone in the parameter space is the topic of our future research.

Figure 7. (a) The quasi-periodic orbit of the swimmer with a single turning cardioid for ${\rm\Delta}{\it\nu}=0.00625$ , and the dispersed chain-like pattern of tracer packets whose advection starts from the tiled elements at $t=0$ . The snapshot of the chain-like chaotic pattern has been taken at $t=625T_{\mathit{slow}}$ . The chain-like pattern is rotating. (b) The swimmer’s periodic orbit with inverted cardioids for ${\rm\Delta}{\it\nu}=0.005$ . Fluid elements are stretched and transported, without chaos, on a curved path along the periphery of the orbit. The snapshots of the four deforming fluid elements have been taken at $t=0$ , $50T_{\mathit{slow}}$ , $100T_{\mathit{slow}}$ , $200T_{\mathit{slow}}$ , and $400T_{\mathit{slow}}$ .