3.1. Swimmer and particle kinematics

We first compare the simulated swimmer locomotion and particle kinematics with results from experiments that have identical set-ups in the swimmer configuration and reservoir size, as described in our previous work (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025). In the experiment, we observed a persistent propulsion over many cycles. As the experimental swimmer was mounted on two vertical rods, we included the rods in one DEM simulation for validation. Figure 2(a) shows the normalised displacement of the swimmer, $y/L$ , over normalised time, $t/T$ . The experimental measurement (dotted curve) and the corresponding simulation result (thick yellow curve) show qualitative agreement, with consistent net locomotion in the $+y$ direction and similar displacement amplitudes per stroke. The simulated net displacement per cycle is larger than that in the experiment. This is likely due to a higher initial packing fraction, $\phi _0=0.60$ , in the simulation, compared with the experiment ( $\phi _0\approx 0.58$ ), as the experimental particle system was air fluidised before each experiment (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025), while particles settle under gravity in DEM. We also plot the result of a simulation without the rods (blue curve), which shows the same $+y$ locomotion direction, a bigger displacement amplitude per stroke, and less net locomotion in each cycle. For simplicity, we will later focus on this simulation in the following sections and refer to it as the base case.

Figure 2.

Swimmer and particle kinematics from the experiment (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025) and from DEM simulations, see configurations in table 1. (a) Swimmer displacement vs time in the experiment (dotted, orange), the DEM result with simulated mounting rods (yellow), the simplified simulation base case (blue), and a frictionless simulation (red). Panels (b) and (c) highlight particles with a displacement larger than five times the average displacement in the DEM results (left panels) and the experimental X-ray CT results (right panels), at $\theta =75^\circ$ during opening and $\theta =25^\circ$ during closing, respectively. Panels (d) and (e) compare the particle displacement field at the mid-plane of the swimmer between DEM results (left panels) and X-ray CT results (right panels), at $\theta =75^\circ$ during opening and $\theta =25^\circ$ during closing, respectively.

Figure 2. Long description

Panel A: A line graph shows swimmer displacement over time. The x-axis represents normalized time (t/T), and the y-axis represents normalized displacement (y/L). The graph includes four data series: experiment (dotted orange line), DEM with rods (yellow line), DEM without rods (blue line), and frictionless DEM (red line). The lines exhibit periodic oscillations, indicating the swimmer’s movement patterns over time. Panel B: Two visuals highlight particles with significant displacement during the opening phase. The left visual shows DEM simulation results, while the right visual shows experimental X-ray CT results. Particles with displacement larger than five times the average are emphasized. Panel C: Two visuals highlight particles with significant displacement during the closing phase. The left visual shows DEM simulation results, while the right visual shows experimental X-ray CT results. Particles with displacement larger than five times the average are emphasized. Panel D: Two visuals compare particle displacement fields at the mid-plane of the swimmer during the opening phase. The left visual shows DEM simulation results, while the right visual shows experimental X-ray CT results. Arrows indicate the direction of particle movement. Panel E: Two visuals compare particle displacement fields at the mid-plane of the swimmer during the closing phase. The left visual shows DEM simulation results, while the right visual shows experimental X-ray CT results. Arrows indicate the direction of particle movement.

To continue our validation, we visualise the particle kinematics using DEM and a quasi-static experiment imaged by X-ray computed tomography (CT) (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025), capturing one full three-dimensional (3-D) scan of the system after an incremental wing rotation. In figure 2(b–e), we first highlight regions with significant particle displacement (panels b and c), and then plot the in-plane ( $x$ – $y$ ) displacement field in the mid-plane of the swimmer (panels d and e), both showing good agreement between the simulation and experiment. Besides the individual particle displacement, we obtained the coarse-grained velocity field of the particles around the swimmer using a binning method, in which the average velocity vector in each bin is weighted by the exact partial volumes of the spheres overlapping with the bin (Strobl et al. Reference Strobl, Formella and Pöschel2016). This allows us to use smaller bins for better spatial resolution, and we further average the binning results over five swimming cycles. Results are displayed in figure 3 for three example angles during the opening and closing strokes, showing the vector field of the $x,y$ component of the particle velocity, along with underlying contours, which represent the strain rate, $\dot {\gamma }$ , calculated as the second invariant of the deviatoric part of the 3-D strain rate tensor.

At intermediate and large opening angles, $\theta =75^{\circ }$ in figure 2(d) and $\theta =50^{\circ }, 71^{\circ }$ in figure 3, each wing induces a rather localised flow in its immediate surroundings, with the radius of the perturbed region comparable to the wing’s length, and it does not overlap with the other wing’s perturbed region. More specifically, the flow around each wing organises into a vortex, in which particles follow the wing rotation as if they are part of the wing. This is similar to the jammed stagnant zones observed in previous studies of penetrating plates (Kang et al. Reference Kang, Feng, Liu and Blumenfeld2018; Feng et al. Reference Feng, Blumenfeld and Liu2019; Harrington et al. Reference Harrington, Xiao and Durian2020), which are cone- or triangular-shaped zones appearing in front of the leading edge. Here, we also see that the strain rate is smaller towards the centre of the wings, reflecting the cone shape. However, we note that, for intruders that are not much larger than the particle diameter, deviatoric deformation and particle rearrangement can still occur (Harrington et al. Reference Harrington, Xiao and Durian2020), and particles do not move as a perfectly rigid cluster. While this adds difficulty in detecting the stagnant zones, we will use force-based methods to better characterise them in § 3.2. With the co-moving stagnant zones, the swimmer-induced flow field is reminiscent of that produced by a counter-rotating vortex pair in a Newtonian liquid, such as the flow generated by a pair of counter-rotating cylinders (van Rees et al. Reference van Rees, Novati and Koumoutsakos2015). While classical fluid dynamics establishes that the propulsion speed of such a vortex pair is linear in the vortex strength (Batchelor Reference Batchelor2000), the analogous behaviour in a granular medium remains unexplored.

Figure 3.

The coarse-grained flow field in the swimmer’s mid-plane, with velocity vectors and shear rate $\dot {\gamma }$ visualised for the frictional base case. The three columns correspond to $\theta = 29^\circ , \, 50^\circ , 71^\circ$ , while the two rows correspond to the opening and closing strokes. The black region covers the bins that overlap with the wings (white). The pink star indicates a highly crowded region and the white star indicates a depleted region.

At smaller angles, $\theta =25^{\circ }$ in figure 2(e) and $\theta =29^{\circ }$ in figure 3, the vortex-like flow field no longer exists, especially during closing, and the two zones perturbed by the wings start to overlap. Near the centre of the swimmer, we see that the in-plane velocity starts to converge, which in turn triggers out-of-plane particle motion, see figure 2(c). In Appendix A, we further show the corresponding velocity fields in three dimensions, confirming a strong upward flow during closing and a downward flow during opening at $\theta =29^\circ$ . At larger $\theta$ , the particle motion is mostly in plane. The flow field at smaller $\theta$ is also more extensive than that at larger $\theta$ , which indicates that localised particle rearrangements become more difficult, signalling intensified jamming effects (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025). We next examine the role of jamming in greater detail using force transmission and structural information.

3.2. Hysteresis in the jammed medium during swimming

In figure 4(a), we visualise the inter-particle contact forces $|\boldsymbol{F}|$ and particle velocities $|\dot {\boldsymbol{r}}|$ in the base case (table 1) at five different $\theta$ during both the opening and closing strokes. For clarity, only the top $0.1\,\%$ of the contact forces are shown, and only particles with velocities $|\dot {\boldsymbol{r}}|\gt 0.33 \omega L$ are shown. During each stroke, stagnant zones with percolated strong forces gradually form near the tips of the rotating wings on the leading edges that press into the medium, as at $\theta =50^\circ$ for both strokes. These contact forces form chains that are rooted from the wings and propagate into the medium, which is a signature of jammed granular media (Bi et al. Reference Bi, Zhang, Chakraborty and Behringer2011; Behringer & Chakraborty Reference Behringer and Chakraborty2018). As each stroke proceeds, a second set of force chains that resist the swimmer’s translation appear near the wing hubs, as seen in $\theta =71^\circ$ for opening and $\theta =29^\circ$ for closing. We refer to these zones as resistive stagnant zones (RSZs) and the other type as propulsive stagnant zones (PSZs), which are separated by the gyration point of a wing that has zero instantaneous velocity in the direction normal to the larger lateral surface of the wings. The distance from a gyration point to a hub is $r_g= |\dot {y}(\theta )|\sin (\theta )/\omega$ , where $\dot {y}$ denotes the swimmer’s translational velocity. As particles within the stagnant zones translate almost rigidly with the swimmer wings, these zones effectively become an extension of the swimming body, and they evolve throughout each stroke. Consequently, hysteresis in the morphology and the strong force density of these stagnant zones between the opening and the closing stroke can break the translational symmetry.

Figure 4.

Force transmission and structural analyses for the base simulation case. (a) Particle velocity and inter-particle contact forces visualised during one swimming cycle. (b) Sketch of spatial decomposition of the swimmer’s surroundings, defining zones I and II. Panels (c) and (d) show the average coordination number $Z$ over $\theta$ for zones I and II, respectively. (e) The $y$ component of the total force exerted on different lateral surfaces of the wings.

To determine the jamming state of the particles near the wings, we decompose the immediate surroundings of the wings into zones I and II, separated by $r_g(\theta )$ , as depicted in figure 4(b). Both zones span the full height of the wings and extend $5d$ from the wings. In each zone, we average the coordination number $Z$ , i.e. the number of contacts per particle, and plot the average as a function of $\theta$ in a swimming cycle in figures 4(c) and 4(d) for zones I and II, respectively. For frictional spheres, the isostatic condition of the coordination number is $(n_{\textit{dim}}+1)\le Z_{\textit{iso}}\le 2n_{\textit{dim}}$ , where $n_{\textit{dim}}=3$ is the dimension (Behringer & Chakraborty Reference Behringer and Chakraborty2018), with the specific value of $Z_{\textit{iso}}$ depending on the friction coefficient and the preparation protocol (Silbert Reference Silbert2010). Throughout the entire cycle, the mean coordination number remains well above the minimal threshold of $n_{ {dim}}+1=4$ . As particles are constantly under the confining pressure from the weight of the particles above them, the particles in zones I and II are always jammed and require finite forces to be rearranged.

In zone I, the coordination number $Z$ increases during opening until $\theta \approx 29^\circ$ , corresponding to the development of force chains in the PSZ in figure 4(a). For $\theta \gt 29^\circ$ , $Z$ plateaus at $Z\approx 6$ as the PSZ becomes fully developed and the force-chain intensity remains unchanged in $50^\circ \leq \theta \leq 80^\circ$ . After the wing-rotation direction changes, $Z$ decreases towards a plateau around $Z=5.7$ , corresponding to the relaxation in the stagnant zones via decompression and particle rearrangement. Similar hysteresis in $Z$ has been reported for jammed granular systems under cyclic shear (Zhang et al. Reference Zhang, Majmudar, Tordesillas and Behringer2010; Bi et al. Reference Bi, Zhang, Chakraborty and Behringer2011; Zhao et al. Reference Zhao, Barés, Zheng, Socolar and Behringer2019, Reference Zhao, Zhao, Wang, Zheng, Chakraborty and Socolar2022).

In zone II, while a similar hysteresis in $Z$ exists, a secondary increase up to $Z\approx 6.5$ occurs near the end of the closing stroke, which is due to the merging of two stagnant zones as in figure 4(a) from $\theta =29^\circ$ to $\theta =20^\circ$ . With $Z\gt 2n_{ {dim}}$ , particles in this region are likely isotropically jammed, whereas particles in zone I are more likely to experience shear jamming, evident by lower $Z$ and anisotropic force chains. During this period, force chains rooted from one wing land on the other, as sketched in $\theta =20^\circ$ , making them effectively ‘internal’ forces of the swimmer with little contribution to the propulsion.

We then separately examine the $y$ component of the forces, $F_y$ , exerted on the four lateral wing surfaces as defined in figure 4(b), which are plotted in figure 4(e) as a function of $\theta$ . While a monotonic increase can be observed in $F_y$ on the back side during opening, the front side $F_y$ plateaus near the end of closing. At this stage, a large $F_y$ arises for the hub as its surrounding region serves as the RSZ. The increasing amount of jammed particles makes it difficult for the swimmer to move in the $-y$ direction, while triggering upward particle flows as in figure 2(c) and Appendix A.

We next vary the DEM configurations to systematically influence the stagnant zone formation and discuss its implication in locomotion, which is realised by simulating different swimmer geometries and a frictionless medium with $\mu =0$ .

For the geometry, we varied the wing gap width $S$ that can influence the intensity of jamming between the wings. To approximate $S=\infty$ , a case with only one wing was simulated. We also altered the swimmer geometry by inserting an inert mid-body to connect the wings to approximate a more realistic-looking swimmer. In figure 5(a), we show how the swimmer’s normalised net displacement per cycle, $\Delta y/L$ , depends on the varied factors. As the gap width increases, $\Delta y/L$ rises until $S/d \approx 6$ , then decreases for $S/d \gt 6$ towards the value for $S/d = \infty$ (the black dotted line). This trend resembles the cooperative effects reported for two closely spaced intruders in granular media, where both the drag force (Agarwal et al. Reference Agarwal, Karsai, Goldman and Kamrin2021a ) and total work (Pravin et al. Reference Pravin, Chang, Han, London, Goldman, Jaeger and Hsieh2021) peak at a separation of approximately $4d$ . Similar to the result with inert rods in figure 2(a), by filling this gap with an inert mid-body, a significant increase in $\Delta y/L$ occurs, accompanied by a stagnant zone (RSZ) with more intense force chains behind the mid-body.

In the simulation with $\mu =0$ , locomotion vanishes (figure 2 a), coinciding with the disappearance of the stagnant zones as seen in the snapshot inside figure 5(a). The initial packing fraction is $\phi _0=0.64$ and coordination number is $Z=6.8$ , which is similar to the simulated value by Silbert et al. (Reference Silbert, Grest and Landry2002), indicating that the frictionless particles are still jammed under gravitational forces. Different from the frictional cases shown in figures 4(a) and 5(a), a larger number of particles are mobilised uniformly around the swimmer, avoiding the hysteretic build-up of strong forces. In other words, the lack of inter-particle friction facilitates local rearrangements that prevent the accumulation of stagnant zones near the wings.

To quantify the ‘intensity’ of stagnant zones induced by the swimmer, we define a strong-contact number, $N_c$ , representing the number of particle contacts where the force magnitude exceeds a threshold $F_c$ . This threshold is set to the top 0.1 % of the contact-force distribution across the entire system, which roughly includes all the visualised strong forces in the stagnant zone in figure 4(a). For the base case with the frictional medium, we extract a threshold of $F_c \approx 1.0$ mN at the time instant $\theta =20^\circ$ , and we apply this value to all other frictional cases. For the frictionless case, we obtained $F_c \approx 0.8$ mN. The threshold value coincides with the largest contact forces at the bottom of the container, due to gravity, as estimated by $\rho \phi _0 g H (\pi d^2/4)\approx 1$ mN, with $H$ being the bed height. In this sense, forces exceeding this threshold are considered ‘strong’ relative to the full granular medium. However, $N_c$ should represent the intensity and localisation of the high-force contact network associated with stagnant zone formation. Thus, this definition requires us to sample near the swimmer, so that forces at the bottom of the container are not included, as they do not contribute to the locomotion of the swimmer. Therefore, we only consider contacts that are within a region, centred around the swimmer, with a size of $55d \times 40d \times 25d$ . The robustness of this definition is examined in Appendix B, where we compare force distributions, normalised high-force tails and various versions of the threshold.

Figure 5.

The influence of jamming on net locomotion. (a) Cycle-averaged swimmer net displacement $\Delta y$ for cases with different gap widths $S$ , with error bars indicating cycle–cycle fluctuations, along with the case with inert mid-body (purple) and the case with $\mu =0$ (black). The black dotted line marks the value of $\Delta y$ for the one-wing case. (b) Cycle-averaged $\Delta y$ vs the contact number difference $\Delta N_c$ . Inset shows the $N_c$ in the base case (black curve with $F_c=1$ mN) and the frictionless case (red curve with $F_c=0.8$ mN) in a swimming cycle. The purple diamond represents the case with inert mid-body, the orange square represents the one-wing case and the green triangle represents the case with $\mu =0$ . The red dashed line is a linear fit, $\Delta y/L = \alpha \ \Delta N_c$ , where $\alpha =8.63\times 10^{-6}$ .

Figure 5. Long description

The image contains two panels, each depicting different types of graphs. Panel A shows a scatter plot with error bars, while Panel B features a scatter plot with a line fit and an inset line graph. Panel A: The scatter plot displays cycle-averaged swimmer net displacement (Δy/L) on the vertical axis against the gap width (S/d) on the horizontal axis. Error bars indicate cycle-to-cycle fluctuations. Various data points are marked with different symbols: blue squares for different gap widths, a purple diamond for the case with an inert mid-body, and a black square for the frictionless case. A black dashed line marks the value of Δy/L for the one-wing case. Panel B: The scatter plot shows cycle-averaged net displacement (Δy/L) on the vertical axis against the contact number difference (ΔNc) on the horizontal axis. Data points are color-coded based on the gap width (S/d) and include a purple diamond for the case with an inert mid-body, an orange square for the one-wing case, and a green triangle for a specific case. A red dashed line represents a linear fit to the data, with the equation Δy/L = 8.63 * 10^(-6) * ΔNc. The inset line graph shows the contact number (Nc) over time (θ in degrees) for the base case (black curve) and the frictionless case (red curve). The black curve has a maximum Nc of 8000, while the red curve has a maximum Nc of 200. The inset also highlights the difference in contact number (ΔNc) between the two cases.

An example of $N_c$ vs $\theta$ in a swimming cycle is shown in the inset of figure 5(b), showing a broken symmetry in the force transmission between the opening and closing strokes, agreeing with the structural signal based on $Z$ in figures 4(c) and 4(d). For the frictionless case, however, a very low number of strong contacts exists near the swimmer, as most of the top 0.1 % forces reside at the bottom, again showing a lack of stagnant zones.

As argued earlier, the excess number of strong contacts towards the end of the closing stroke can limit the swimmer’s translation. Therefore, we quantify the degree of stagnant-zone-induced symmetry breaking by calculating the difference of $N_c$ between the fully closed and fully open states, yielding

(3.1) \begin{align} \Delta N_c = N_c(\theta _{\mathrm{c}})-N_c(\theta _{\mathrm{o}}). \end{align}

Figure 5(b) shows the net displacement $\Delta y/L$ with respect to $\Delta N_c$ , revealing a linear relation between these two physical quantities across all cases with varied swimmer geometry and friction. This quantitatively demonstrates how stagnant-zone-induced effects in jammed media break the symmetry between the two reciprocal strokes of the scallop swimmer and enable net locomotion in the $+y$ direction.

To rationalise the relation between $\Delta N_c$ and $\Delta y/L$ , we adopt an argument based on ‘kinematic compatibility’ (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025). For each stroke, the swimmer’s wings sweep across a volume that is $L^3(\theta _{{o}}-\theta _{{c}})$ . If the particles rearrange freely, this volume can be compensated by the swimmer translating a distance $y_s$ accordingly, sweeping a volume that is approximately $2L^2\sin (\theta _{{m}})y_s$ , with $\theta _{{m}}=(\theta _{{o}}+\theta _{{c}})/2$ . Equating the two volumes, we obtain $y_s\approx 0.7L$ , which is comparable to the observed $y_s$ in figure 2(a), especially for the frictionless case that has the largest $y_s$ per stroke. When particles can easily rearrange, e.g. with $\theta =50^\circ$ in figure 3, the vortex-like flow field develops, which corresponds to the flow of particles from zone I towards zone II (figure 4 b) during opening and vice versa.

For swimming in a jammed, solid-like medium with frictional particles, a subset of particles, characterised by $N_c$ , develop strong resistance and do not rearrange to accommodate the swimmer’s translational motion. This effect is most pronounced at small $\theta$ during the closing stroke when $N_c$ peaks. For example, at $\theta =29^\circ$ , the vortex-like particle transport breaks down on the inner side of the wings. Particles in zone II (pink star in figure 3) are highly crowded, as indicated by the converging in-plane velocity and the concentration of strong forces (figure 4 a), and they cannot effectively flow towards zone I (white star). Without these particles replenishing zone I, the hubs moving away from zone I (in the $-y$ direction) would tend to create a ‘vacuum’ in the zone. This is kinematically incompatible in a dense granular medium, and therefore the swimmer translation $y_s$ is hampered. Instead of pumping particles in the $y$ direction to propel itself, the swimmer lifts particles vertically in the $z$ direction, as shown in Appendix A.

A simplified way to account for the effect of $N_c$ is to consider the corresponding characteristic volume, $2N_cd^3/\bar {Z}$ , as a ‘negative’ volume to be subtracted from the wing-rotation-swept volume that is available for the swimmer translation. Here, $\bar {Z}$ is a nominal coordination number, introduced to convert the contact count to the particle count. This viewpoint is equivalent to the viewpoint that the swimmer is effectively enlarged by $2N_cd^3/\bar {Z}$ and needs to move together with these stagnant particles. The characteristic volume balance for each stroke is then $L^3-2N_cd^3/\bar {Z}=2L^2y_s$ , noting that the $\theta$ -related factors are dropped and the $N_c$ used here is measured at the end of each stroke, consistent with the calculation of $\Delta N_c$ . Taking the difference between the two strokes and estimating $\bar {Z}=6$ , we recover the linear relation, $\Delta y/L=[d^3/(6L^3)]\Delta N_c$ , noting that the rotation-swept volume cancels out due to reciprocity. Here, the proportionality factor, $d^3/(6L^3)=5\times 10^{-5}$ , is within an order of magnitude of the fitted constant, $\alpha =8.63\times 10^{-6}$ , from figure 5(b), which is reasonable given the simplicity of the argument. For instance, not all $N_c$ may contribute to such ‘negative’ volume, and accounting for that would lower the proportionality factor. That said, this linear relation may break down at shallower burial depths, where particles are less confined and the free surface is free to deform plastically, altering the history-dependent response of the medium.

3.3. Role of swimmer inertia in reciprocal swimming

A second contribution to the locomotion, associated with the swimmer’s inertia, emerges with increasing swimmer mass $m_s$ and decreasing swimming cycle period $T$ . To quantify the inertial effect, we consider a coasting time

(3.2) \begin{align} T_{c} = \frac {v}{F_D/m_s} = \frac {L/T} {P_sL^2/m_s}=\frac {m_s} {P_sLT}, \end{align}

where $L/T$ sets a characteristic swimmer translation speed that is comparable to the swimmer’s speed near the end of a stroke, i.e. the slope of the curves in figure 2(a). We consider a characteristic acceleration, $P_sL^2/m_s$ , where $P_s=\rho \phi _0 gH_0$ is a hydrostatic-like pressure, and $H_0$ is the swimmer depth with respect to the free surface. In this case, $P_sL^2$ serves as a normalisation factor, which should be of the same order of magnitude as the drag force. This term comes from the rate-independent term in a previously established drag law in granular intrusion, $F_D = a_1 P_s A + a_2\rho \phi _0 A v^2$ , where $A\sim L^2$ is the projected frontal area, $v\sim L/T$ is the intruder speed and $a_1, a_2$ are material and geometric prefactors (Albert et al. Reference Albert, Pfeifer, Barabási and Schiffer1999; Katsuragi & Durian Reference Katsuragi and Durian2007; Clark et al. Reference Clark, Kondic and Behringer2012; Brzinski III et al. Reference Brzinski III, Mayor and Durian2013; Hilton & Tordesillas Reference Hilton and Tordesillas2013; Takehara & Okumura Reference Takehara and Okumura2014). The Froude number ${Fr}= v/\sqrt {gH_0} = {v}/{\sqrt {P_s/(\rho \phi _0)}} \lesssim 0.54$ across the entire parameter sweep of table 1, so the rate-independent term dominates and $F_D\approx a_1 P_s L^2$ . Thus, the use of $P_sL^2$ as a normalisation factor is reasonable and convenient, as it remains a constant across all our simulation configurations. Under this definition, $T_c$ reflects a characteristic time for a moving swimmer to ‘coast down’ until its translational motion ceases due to the medium’s drag (Hubert et al. Reference Hubert, Trosman, Collard, Sukhov, Harting, Vandewalle and Smith2021), were it to stop actuating.

Figure 6.

Locomotion in the dynamic regime. (a) The normalised swimmer displacement $y$ vs time in the quasi-static (dashed blue curve) and dynamic (solid green curve) regimes. The quasi-static case is with $m_s=0.60 \text{ g}$ , $T=2.0 \text{ s}$ and $T_c/T=5.16\times 10^{-5}$ , and the dynamic case has $m_s=2.38 \text{ g}$ , $T=0.05 \text{ s}$ and $T_c/T=0.33$ . (b) Normalised $y$ vs $t/T$ for one swimming cycle for the quasi-static (dashed curve) and dynamic (solid curve) cases. The grey region marks the opening stroke, and the vertical dashed lines mark the time instants at which the direction of the swimmer’s translational motion switches. (c) Normalised $y$ vs $\theta$ for the two cases, with lighter colours representing opening and darker colours representing closing. The horizontal dashed green lines mark the $y$ at which the direction of the swimmer’s translational motion switches, while the dashed grey lines mark the direction switch of the wing rotation.

Figure 6(a) shows the swimmer displacement over time for a simulation in the dynamic regime with $T_c/T=0.33$ , along with the base case with $T_c/T=5.16\times 10^{-5}$ . Contrary to the forward ( $+y$ ) quasi-static locomotion (Xiao et al. Reference Xiao, Torres, Sack and Pöschel2025), a persistent locomotion occurs in the $-y$ direction, indicating a distinct mechanism in this regime. This direction agrees with the theoretical prediction of a scallop-like swimmer with non-negligible inertia in Newtonian fluids at low Reynolds numbers (Gonzalez-Rodriguez & Lauga Reference Gonzalez-Rodriguez and Lauga2009), despite the differences in the nature of forcing and problem set-up.

We scrutinise $y/L$ within a single swimming cycle, as shown in figure 6(b). As the swimmer never stops actuating, there is no strict coasting. Instead, we see a retardation process characterised by a time interval, $\Delta t^i_R$ , measuring the delay between the instant when the wing rotation reverses and the instant when the swimmer’s translation direction reverses (marked by the vertical dashed lines in figure 6 b), which can be seen after both the opening ( $i=\mathrm{o}$ ) and closing ( $i=\mathrm{c}$ ) strokes. Likewise, we measured the corresponding retardation distance, $\Delta y^i_R$ , which is made clear in figure 6(c) that shows $y$ against the wing angle $\theta$ in a cycle. Both results show that, for $T_c/T=0.33$ , significant retardation occurs, confirming that the case is in the dynamic regime. For $T_c/T=5.16\times 10^{-5}$ , we see that both $\Delta t^i_R$ and the corresponding $\Delta y^i_R$ are negligible, indicating a quasi-static regime where the swimmer inertia is insignificant.

For $T_c/T=0.33$ , the retardation after the closing stroke is more prolonged in comparison with the retardation after the opening stroke, yielding $\Delta t^c_R\gt \Delta t^o_R$ and $|\Delta y^c_R|\gt |\Delta y^o_R|$ . To further understand this difference, a parametric study was performed, varying $T$ and $m_s$ to determine the dependence of $\Delta y^i_R$ and $\Delta t^i_R$ on different $T_c/T$ in both frictional and frictionless media. The absence of friction eliminates the formation of stagnant zones, as discussed in § 3.2, which helps highlight the role of inertial coasting. We show $|\Delta y^c_R|-|\Delta y^o_R|$ in figure 7(a), and $\Delta t^c_R-\Delta t^o_R$ in figure 7(b), both confirming that the retardation is more prolonged after the end of the closing stroke, and the differences increase with $T_{c}/T$ .

Figure 7.

The dependence of $\Delta y^i_R$ and $\Delta t^i_R$ on different $T_c/T$ in both frictional and frictionless media by varying $T$ and $m_s$ according to table 1. (a) The difference in retardation distance between the two strokes vs the coasting time $T_c/T$ in various cases with frictional (blue circle) and frictionless (red triangle) media. Panel (b) is the same as (a) but for the difference in the retardation time. (c) The contact number difference $\Delta N_c$ vs the coasting time $T_c/T$ for the cases in the frictional medium. The two snapshots of the swimmer with force chains visualised correspond to the two cases with $T_c/T=0.005$ and $T_c/T=0.33$ , which show the intensity of the contact-force network in a similar colour scheme as figure 4(a).

Figure 7. Long description

Three scatter plots depict the dependence of various parameters on coasting time in frictional and frictionless media. Panel A: A scatter plot shows the difference in retardation distance between two strokes versus coasting time. The x-axis represents the coasting time (T_c/T) and the y-axis represents the difference in retardation distance (|Δy_R^c| - |Δy_R^o|)/L. Data points are marked with blue circles for frictional media and red triangles for frictionless media. Panel B: Another scatter plot shows the difference in retardation time versus coasting time. The x-axis represents the coasting time (T_c/T) and the y-axis represents the difference in retardation time (Δt_R^c - Δt_R^o)/T. Data points are marked similarly with blue circles for frictional media and red triangles for frictionless media. Panel C: A scatter plot shows the contact number difference versus coasting time for cases in the frictional medium. The x-axis represents the coasting time (T_c/T) and the y-axis represents the contact number difference (ΔN_c). Blue circles represent data points, and two snapshots of the swimmer with force chains are visualized, corresponding to two cases with different coasting times. The intensity of the contact-force network is shown in a similar color scheme as figure 4(a).

The difference in the retardation behaviour is reminiscent of the result in Hubert et al. (Reference Hubert, Trosman, Collard, Sukhov, Harting, Vandewalle and Smith2021) that a floating dumbbell swimmer can propel itself provided that the respective coasting times for its reciprocal strokes are different. Likewise, the symmetry-breaking mechanism in the inertial regime of the scallop swimmer is related to this difference in the retardation behaviour. In the reciprocal strokes of opening and closing, the swimmer gains similar peak velocity, as reflected by the slope of $y$ in figure 6(c) near $\theta =50^\circ$ , noting that $\theta$ is linear with $t$ . However, near the end of each stroke, including the beginning period of the subsequent stroke, the resistive drag force that opposes the swimmer’s translation is different between opening and closing. This force should be proportional to the swimmer’s exposed area in the $x$ – $z$ plane, $2L^2\sin (\theta )$ , which is larger after opening than closing. This results in the more prolonged retardation after the closing stroke, hence the preferred backward ( $-y$ ) locomotion direction.

Figure 8.

Particle-scale analyses of the frictionless granular medium at two different $T_c$ . (a) Visualised particle velocity and inter-particle contact force at the end of each stroke and at $0.02 \,T$ into the subsequent stroke. The first row is for a swimmer with $m_s=0.60 \text{ g}$ , $T=2.0 \text{ s}$ and $T_c/T=5.16\times 10^{-5}$ , and the second row is for $m_s=2.38 \text{ g}$ , $T=0.05 \text{ s}$ and $T_c/T=0.33$ . Only the top 0.1 % of the contact forces and the particles with $|\boldsymbol{\dot r}|\gt 0.33\omega L$ are shown. (b) Swimmer displacement over time. (c) Mass of mobilised particles, whose velocity is larger than $0.33 \omega L$ during a swimming cycle. (d) Momentum magnitude of the mobilised particles. (e) The $y$ component of the medium’s resistive force applied to the swimmer. In (b–e), the red dashed curves correspond to $T_c/T=5.16\times 10^{-5}$ and the green solid ones to $T_c/T=0.33$ . The lighter and darker curves correspond to the opening and closing strokes, respectively.

To isolate the contribution of inertia to the backward locomotion, we further examine particle-scale features in the frictionless medium, which suppresses the friction-induced hysteresis as in § 3.2. We compare a quasi-static case with $T_c/T = 5.16\times 10^{-5}$ and a dynamic one with $T_c/T = 0.33$ , visualising the inter-particle contact forces and particle velocities at the end of each stroke and at $0.02 \,T$ into the subsequent stroke, which is within the retardation process for the dynamic case. Without particle friction, the backward locomotion displacement over time is larger for $T_c/T = 0.33$ than that with frictional particles, which can be seen by comparing figure 8(b) with figure 6(a). We also track the total mass, $M_{{mob}}$ , of mobilised particles (with $|\boldsymbol{\dot r}|\gt 0.33\omega L$ ), the corresponding total particle momentum magnitude $P_{\textit{mob}}(t) = \sum _j m_j|\dot {\boldsymbol{r}}_j|$ and the $y$ component of total resistive force on the swimmer, which are shown in figure 8(c–e), respectively.

Between the two cases, major differences are found in the particle-scale response to the wing-rotation reversal, as shown in figure 8(a). For the quasi-static case, almost no concentrated, strong-contact network forms at $t/T=\{0.0, 0.5\}$ and at $t/T=\{0.02,0.52\}$ . A wide range of particles are in motion and the velocity of the mobilised particles is smaller at $t/T=\{0.02,0.52\}$ compared with those at $\{0.0, 0.5\}$ . As $T_c/T$ is small, particles have enough time to reorganise and dissipate the transferred momentum from the swimmer, as confirmed in figure 8(c), showing a persistently large mobilised mass over a cycle. Also, the transferred momentum and resistive force in this regime are negligible, as shown in figure 8(d,e), confirming that no granular inertia is carried across strokes.

For the dynamic case, however, a strong, concentrated contact network emerges after $t/T=\{0.0, 0.5\}$ and further intensifies at $t/T=\{0.02,0.52\}$ . This is a consequence of the swimmer coasting, i.e. the wings both rotating into and translating into the corresponding particles at high speeds. As a result, particles do not have enough time to rearrange to dissipate the transferred momentum throughout the medium and instead develop large contact forces. In fact, far fewer particles exceed the $0.33\omega L$ speed threshold, as shown qualitatively in figure 8(a) and quantitatively in figure 8(c). Similarly, we observe an increase in momentum and resistive force at $\theta$ that corresponds to $t/T=\{0.02,0.52\}$ in figure 8(d,e). As a result, all three quantities exhibit peaks during swimmer coasting. More importantly, these peaks are asymmetric between the two strokes, with peaks following the opening stroke being larger in magnitude compared with those following the closing stroke, which, again, is due to the geometric asymmetry of the wings. At $\theta \approx 20^\circ$ , the concentrated strong forces result in mostly lateral forces (in the $x$ direction) on the wings, whereas at $\theta \approx 80^\circ$ , $|F_y|$ is larger, which impairs the coasting motion as discussed earlier.

Even without friction, particles with concentrated strong forces (figure 8 a) also constitute stagnant zones due to the lack of relaxation, i.e. their jammed packing structure remains ‘stagnant’. However, key differences exist between these zones that are rate-dependent and the friction-dependent zones, as discussed in § 3.2. The rate-dependent zones appear at the beginning of each stroke and resist coasting motion in the direction set by the previous stroke, whereas the friction-dependent zones become most prominent at the end of each stroke and resist motion in the direction set by the current stroke. Hence, the respective propulsion mechanisms, related to the symmetry breaking of these zones, result in opposite propulsion directions. Of course, for dynamic swimming in a frictional medium (figure 6), both mechanisms coexist, and we next seek a more unified description.

3.4. Coexistence and unified scaling of the two mechanisms

Having characterised the jamming mechanism in § 3.2 and the inertial mechanism in § 3.3 in isolation, we now examine how the two mechanisms coexist based on the simulations of figure 7. The increase of $|\Delta y_R^c|-|\Delta y_R^o|$ and $\Delta t_R^c-\Delta t_R^o$ with $T_c$ is linear in the frictionless cases, while it saturates in the frictional cases at large $T_c$ . The saturation is likely due to the competition from the stagnant zone formation due to jamming, as discussed in § 3.2. For example, during retardation after the closing stroke, the rotating wing is pressing towards the $-y$ direction while the swimmer is also moving towards the $-y$ direction under its inertia. This may significantly compress the particles at the back side of the swimmer, inducing a large RSZ that hinders the retardation distance and results in the saturation in figures 7(a) and 7(b). From example snapshots inside figure 7(c), we see a larger RSZ for the case with higher $T_c/T$ .

Therefore, the jamming-based and inertia-based mechanisms are not mutually exclusive regimes; they coexist and contribute to locomotion with different parameter dependencies. The jamming-based contribution is governed by the friction coefficient $\mu$ and the hysteretic nature of the medium, which manifests as asymmetric stagnant zones towards the end of the two strokes, as captured by $\Delta N_c$ . The inertial contribution is governed by the ratio of the coasting time to the swimming period, $T_c/T$ , which is more relevant to the swimmer dynamics at the beginning of the two strokes. It is also worth noting that the two mechanisms mentioned in this study appear as competing contributions in the net locomotion at high $T_c/T$ . While the jamming-based mechanism favours the forward motion, the inertia-based mechanism favours the backward motion.

We then seek to construct a unified description considering both mechanisms. In figure 9(a), we relate the net swimmer locomotion to the coasting time $T_c/T$ , which reflects both the significance of the swimmer’s inertia and the difference between the coasting after the closing and opening strokes. In the frictional cases, a plateau exists at low $T_c/T$ , showing a quasi-static regime as discussed in § 3.2. At $T_c/T\approx 4 \times 10^{-3}$ , a transition to a dynamic regime occurs, showing a decrease in forward locomotion or even the presence of backward locomotion. In the frictionless cases, the same plateau exists in the quasi-static regime, but without forward locomotion. A transition occurs at a similar $T_c/T$ as in the frictional cases, after which the swimmer shows a significant increase in backward locomotion. The difference in $\Delta y/L$ between the frictional and the frictionless cases is related to $\Delta N_c$ , as shown in figure 5(b). After subtracting $\alpha \Delta N_c$ from $\Delta y/L$ , the frictional and frictionless results collapse in figure 9(b), where a master power-law curve fits well to all data, as

(3.3) \begin{align} \Delta y / L = \alpha \Delta N_c - \beta \ (T_c/T)^{n} , \end{align}

with the fitted $\beta =0.24$ and $n=0.44$ , highlighting the respective roles of jamming and swimmer inertia in locomotion. The collapse of the data suggests that the locomotion of the scallop swimmer can be decomposed into a contribution from jammed stagnant zone formation, characterised by $N_c$ , and a contribution from the inertial coasting, characterised by $T_c/T$ . We note that (3.3) should not be extrapolated to $T_c/T=\infty$ . As the swimmer keeps increasing its flapping frequency, the granular medium can be strongly agitated and transition to a gaseous behaviour, which can possibly lead to a decrease in $|\Delta y|$ . Similarly, if the swimmer’s mass $m_s$ keeps increasing, the swimmer may not be able to gain enough momentum from a stroke, also resulting in decreased $|\Delta y|$ .

Figure 9.

Quantified relations between swimmer locomotion and inertial coasting. (a) Normalised swimmer displacement vs the coasting time in the frictional medium (blue circles) and the frictionless medium (red triangles). Green circle markers show the simulation cases demonstrated in figure 6 and the green triangle marker shows the case in the frictionless medium demonstrated in figure 2(a). The vertical dotted line indicates the transition from the quasi-static regime to the dynamic regime. (b) The same plot but with $\Delta y/L$ offset by the jamming contribution $\alpha \Delta N_c$ . The dashed curve is the fitted (3.3).