Snapping elastic disks as microswimmers: swimming at low Reynolds numbers by shape hysteresis

Abstract

We illustrate a concept for shape-changing microswimmers, which exploits the hysteresis of a shape transition of an elastic object, by an elastic disk undergoing cyclic localized swelling. Driving the control parameter of a hysteretic shape transition in a completely time-reversible manner gives rise to a non-time-reversible shape sequence and a net swimming motion if the elastic object is immersed into a viscous fluid. We prove this concept with a microswimmer which is a flat circular elastic disk that undergoes a transition into a dome-like shape by localized swelling of an inner disk. The control parameter of this shape transition is a scalar swelling factor of the disk material. With a fixed outer frame with an additional attractive interaction in the central region, the shape transition between flat and dome-like shape becomes hysteretic and resembles a hysteretic opening and closing of a scallop. Employing Stokesian dynamics simulations of a discretized version of the disk we show that the swimmer is effectively moving into the direction of the opening of the dome in a viscous fluid if the swelling parameter is changed in a time-reversible manner. The swimming mechanism can be qualitatively reproduced by a simple 9-bead model.

---

1 Introduction

Swimming on the microscale at low Reynolds numbers requires special propulsion mechanisms which are effective in the presence of dominating viscous forces. Phoretic swimmers create gradients in external fields such as concentration or temperature which in turn give rise to symmetry-breaking interfacial forces leading to propulsion if they overcome the friction force of the microswimmer.¹ Besides phoretic swimmers, self-deforming or shape-changing swimmers are the largest class of microswimmers. They deform their body in a cyclic way in order to propel. This general principle has the advantage that it works independently of the environment, i.e., it does not require an external field, that eventually changes the properties of the fluid. The disadvantage of swimming by deformation is, on the other hand, that there are necessarily “moving parts” causing additional viscous flow in the fluid and forces on the swimmer. As a consequence, at low Reynolds numbers, the cyclic deformation pattern must not be invariant under time-reversal: the scallop theorem formulated by Purcell states that periodic reciprocal patterns of deformation can not lead to an effective net motion on the microscale because of the linearity of the Stokes equation.²

In nature, many different examples of deformation swimmers can be found such as bacteria, algae and spermatozoa.^3,4 These natural swimmers often rely on the movement of a few flagella or many cilia on their surface.^5–8 Flagella employ a periodic forcing but overcome the scallop theorem by exploiting friction along the elastic flagellum to break time-reversibility. This requires a matching of driving and frictional damping time scales for efficient propulsion. Often it also requires the ability of local actuation for the periodic forcing. This makes this concept hard to reproduce or imitate in a controlled fashion in an artificial system.^9,10

Another basic strategy to overcome the scallop theorem are deformation cycles that involve at least two control parameters and drive the swimmer periodically along a closed contour in this at least two-dimensional parameter space. Different shape changing artificial swimmers have been developed based on this concept starting with Purcell's three-link swimmer² and including swimmers performing small deformations of spheres, circles, or cylinders,^11–15 or shape-changing vesicles.¹⁶ The most simple shape changing microswimmer is arguably the one-dimensional linear three-bead swimmer developed by Najafi and Golestanian, where three beads change their distance in a non-time-reversible way.^17,18 By extending the linear three bead arrangement to a second dimension in a triangular shape, a three bead swimmer can perform two-dimensional motions (circles)¹⁹ and steer.²⁰ Nevertheless, despite the simplicity of the concept, this type of swimmer is difficult to implement experimentally because it requires fine control over, at least, two control parameters such as the bead positions of the three-bead swimmer.^21,22

We employ a different general strategy in order to overcome the scallop theorem, which is widely applicable and only involves control of a single global and scalar control parameter, which couples, however, to a hysteretic (or bistable) shape transition of the system, see Fig. 1. If also the sequence of shapes exhibits hysteresis, this converts the time-reversible motion in one-dimensional control parameter space into a non-time-reversible motion in a higher-dimensional parameter or shape space. Hysteretic shape transitions can be realized, for example, by using the intrinsic properties of elastic materials. In this work, we will realize such a hysteretic shape transition based on a swelling process of a flat and thin circular elastic disk, where material swelling with swelling ratios of only a few percent in the central region of the disk leads to a shape transition from the flat disk shape into curved conformations, such as a dome-like shape.^23–25 The snapping into an elliptic dome-like shape actually faintly resembles the opening and closing of a scallop. By further enhancing the elastic disk with a fixed frame with attractive interactions we can endow this transition with genuine hysteretic effects. These hysteresis effects allow us to break the reciprocity of the shape cycle although we employ simple cyclic and fully time-reversible oscillations of the swelling factor as single global and scalar control parameter. The main point of this paper is to give the proof of concept that this leads to net propulsion and is a viable realization of a microswimmer. The principle of exploiting a periodically driven hysteretic shape transition in order to achieve net propulsion has been introduced in ref. 26 using the buckling transition of spherical elastic shells as propulsion mechanism. The buckling of spherical shells is a subcritical hysteretic shape transition,^27–29 which will turn out to be conceptually similar to the elastic instability triggered in the elastic disk by localized swelling.

Fig. 1 A completely time-reversible oscillation of the control parameter (horizontal axis, here: swelling ratio) gives rise to non-time-reversible shape cycle of an elastic disk because of hysteresis of the triggered shape transition. The deformation cycle between shapes A and C resembles the opening and closing of a scallop but is hysteretic (compressed shape B and the transition state between B and C are only visited upon swelling).

Deformation cycles of elastic materials have been applied before to design artificial swimmers. In ref. 30, structured light fields were used to drive elastic deformation waves on swimmers with a homogeneous body made of a soft material. Other approaches focused on elastic double layers, where swelling of one layer can induce bending; such externally controllable swelling layers can be engineered, for example, using thermoresponsive microgels.³¹ These ideas were used to design swimmers with a helical structure that can be propelled by conformation changes of the helix.^31–33 Here conformation changes are non-reciprocal, partly because of hysteresis effects in the heating cycle of the thermoresponsive gel.

Deformation swimmers are relatively slow in general, because the swimming distance only scales quadratically with the deformation displacement for many deformation swimmers.^11,12 Therefore, a high frequency of conformation changes is needed in order to achieve a significant swimming velocity. This applies also to the concept of swimming by hysteretic shape changes. The driving frequency is, however, not limited by an additional damping time scale (as, for example, in flagellar motion) because breaking of the time-reversibility is inherent in the hysteretic shape sequence itself but, at high driving frequencies, one could leave the realm of low Reynolds numbers.

The paper is organized as follows. At first, we present a “dry” numerical analysis of the elastic deformation cycle of the elastic disk with material swelling in the interior in the absence of a surrounding fluid and corresponding hydrodynamic interactions. We quantify the hysteretic effects of the swelling transition and show that there is most likely no genuine hysteresis in experiments on simple swelling disks. We can enhance the model system by adding an additional fixed frame with an attractive interaction to the disk in order to generate a robust and genuine hysteretic shape transition. Then we perform a “wet” hydrodynamic simulation featuring a Rotne–Prager interaction in order to proof the swimming ability of this system and characterize the net propulsion. Finally, a simplified 9-bead model is presented that mimics the essence of the underlying swimming mechanism and is able to qualitatively reproduce and explain its main characteristics.

2 Swelling of a flat disk

2.1 Theory

In the following, we consider a flat circular disk with radius R_out and thickness h. The disk shall be very thin (h ≪ R_out), so we can use a two-dimensional model. We parametrize our two-dimensional disk in polar coordinates (r,φ). The basic idea is to deform the disk into a curved shape by a localized, i.e., inhomogeneous swelling process which changes the disk's metric.^23–25 By swelling we mean a local isotropic swelling, where the rest lengths of fibers change by a position-dependent factor A(r,φ) independent of fiber orientation. In the following, we restrict ourselves to radially symmetric swelling functions A(r); the neutral case of the flat disk is represented by A(r) = 1, A(r) > 1 corresponds to local swelling, A(r) < 1 to local shrinking of fibers. In order to calculate the change in metric by a swelling function A(r) we re-parametrize the deformed shape using Gaussian normal coordinates. The Gaussian radial coordinate is given by , which is the distance of a point to the origin of the coordinate system following the surface and reduces to the standard radial coordinate r for a flat disk A(r) = 1. The angular coordinate φ remains unchanged. In Gaussian coordinates a deformed fiber in radial direction has length dl_ρ = dρ = A(r)dr, a deformed circumferential fiber a length dl_φ = rA(r)dφ. In general, the fiber length is related to the metric by dl² = g_ρρdρ² + 2g_ρφdρdφ + g_φφdφ² such that the metric tensor of the swollen disk can be read off asin Gaussian normal coordinates. This is the so-called target metric which represents the preferred equilibrium state of the swollen disk.²⁴ According to the Theorema Egregrium the Gaussian curvature K̄ can be deduced solely from the metric tensor. Using the Brioschi formula (with respect to the Gaussian coordinates ρ and φ of the metric) we findas a function of the Gaussian radial coordinate ρ. We can transform to the standard radial coordinate r by using dr/dρ = 1/A(r), which gives

This is the Gaussian curvature if a shape with the metric (1) could be embedded into three-dimensional Euclidian space.

In order to deform the disk into a curved shape, we define a simple class of swelling patterns. Similar to Pezulla et al.,²⁵ we divide the disk into two parts: An inner disk with 0 ≤ r ≤ R_in and an outer annulus with R_in < r ≤ R_out. Within these two regions, the swelling function A(r) shall be piecewise constant. To simplify things further, we define the inner disk to swell with a constant factor α, while the outer annulus shall always do the exact opposite, i.e., it shrinks with the inverse constant factor 1/α. In total, the considered swelling functions A(r) can be written as

The swelling process is thus defined by two simple control parameters: the swelling factor α and the geometrical ratio R_in/R_out, which will be kept constant at R_in/R_out = 0.5 in the following so that we can focus on the influence of α. We can distinguish between two general cases: (a) α > 1 means that the material in the inner disk is swelling, while the outer annulus is shrinking, which is illustrated in Fig. 2(a). The opposite case, (b) α < 1, leads to a shrinking inner disk and a swelling annulus (Fig. 2(b)). In order to apply the Brioschi formula (3), we have to smear out the step function at r = R_in in eqn (4). For α > 1, we have A′(r) = 0 except around r = R_in, where −A′(r) is peaked. Therefore, the last term in the denominator in eqn (3) dominates and gives a positive Gaussian curvature which is peaked around r = R_in if a surface with the piecewise metric (1) could be embedded into three-dimensional space. Although the metric (1) is piecewise flat, Gaussian curvature has to be introduced because the inner disk is bonded to the outer annulus. This Gaussian curvature is positive for α > 1 and negative for α < 1. Because the metric can actually not be embedded, the Gaussian curvature will be redistributed on the entire disk to minimize the elastic energy. We thus expect a target elliptic shape with an overall positive Gaussian curvature for α > 1, and a hyperbolic target shape with a negative Gaussian curvature for α < 1.

Fig. 2 Schematic demonstration of the swelling pattern. The disk is divided into two parts. The outer annulus (blue) swells with a constant factor α, the dinner disk (red) with 1/α. (a): With α > 1, the inner disk expands and the annulus shrinks. (b): Shrinking in the inner region and material expansion in the annulus with α < 1.

Obviously, for α ≠ 1 the outer annulus and the inner disk are incompatible at r = R_in. Therefore, the surface described by the metric (1) has no immersion in the Euclidean three-dimensional embedding space and is an example of non-Euclidean geometry.²⁴ The actual shape of the surface in three-dimensional space is then defined by the minimization of the elastic energy, where the elastic energy of the deformed swollen state is defined with respect to the above target metric ḡ from eqn (1). The incompatibility of the two parts and the resulting non-existence of an immersion means that there always is a residual elastic energy after minimization.²⁴

The elastic energy contains a stretching and a bending contribution. The stretching contribution is caused by strains ε_ij = (g_ij − ḡ_ij)/2, where g is the actual metric that the deformed swollen state assumes, and given by²⁴

where we use Einstein summation, raising of indices is performed with the target metric, and (u₁,u₂) = (ρ,φ) in Gaussian normal parametrization. The elastic tensor A^ijkl is given by the two dimensional Young modulus Y_2D and the Poisson ratio ν, which characterize the stretching elasticity of the disk material. The stretching energy thus penalizes deviations from the target metric.

Likewise, the bending energy is defined with the curvature tensor L and with respect to a target curvature tensor L̄, which represents a spontaneous curvature of the material. We assume that local isotropic swelling does not introduce any spontaneous curvature to the system such that L̄ = 0. The general expression for the bending energy is²⁴

where the last line applies to L̄ = 0, H is the mean curvature, K the Gaussian curvature, and κ_B = Y_2Dh²/(12(1 − ν²)) the bending modulus of the disk. The bending energy penalizes deviations from the flat shape for vanishing target curvature L̄ = 0. The last line in (8) is an approximation because we assume ḡ_ij ≈ g_ij in A^ijkl. Typical strains ε_ij are ∝ (1 − α)² such that corrections are O((1 − α)²E_B) and will be small at the transition for thin disks (see eqn (17) and Fig. 6 below).

For numerical energy minimization the disk and its elastic energies (5) and (7) have to be suitably discretized.

2.1.1 Model. We calculate the disk's shape with the help of a numerical energy minimization and use a simple spring mesh model for discretization. The disk is triangulated with a Delaunay triangulation (implemented with the fade2D library³⁴), where every edge i between two vertices represents a mechanical spring with a rest length l_i. The fineness of the mesh is controlled by the number of vertices n_B on the boundary of the disk. In this model, a swelling process is performed by a simple multiplication of the springs' rest lengths with the swelling function A(r). The discretized version of the elastic stretching energy (5) can be written as the sum over all spring energies,

The vectors r⃑_2,i and r⃑_1,i describe the positions of the vertices that define the beginning and the end of a spring. The spring constants are denoted by k_i. In a hexagonal mesh, the two-dimensional Young modulus Y_2D is given by the spring constant k and the Poisson ratio ν is fixed,^35,36

In order to evaluate the bending energy (8) on the spring mesh, the curvatures H and K have to be calculated on the mesh. The mean curvature H_i at a mesh vertex i can be expressed in terms of an area gradient:³⁷

The quantity A_i represents the area in the mesh that is associated to the vertex i, see the colored area in Fig. 3. The gradient ∇_iA_i then describes derivatives of this area with respect to the coordinates of the vertex i. The Gaussian curvature K, on the other hand, can be calculated using the Gauss–Bonnet-theorem. We find

where θ_j is the angle between the neighboring vertices j and j + 1 of the vertex i located at r⃑_i, see Fig. 3. Finally, the discretized bending energy (8) becomes

Fig. 3 Illustration of the direct neighborhood of a vertex at position r⃑_i. The area A_i that is associated with this vertex is shown cyan. The angle between the springs to the neighbor vertices at r⃑_j and r⃑_j+1 is called θ_j and is used in the calculation of the Gaussian curvature.

The total energy E = E_s + E_B has to be minimized with respect to all vertex coordinates in the three-dimensional embedding space. In order to overcome possible local energy minima, small fluctuations can be added to the vertex coordinates in terms of a random displacement with . After minimizing the global energy minimum with respect to all vertex positions, the resulting mesh represents the preferred configuration of the swollen and deformed disk, see the illustration in Fig. 4.

Fig. 4 Example Delaunay triangulation with n_B = 40 vertices on the boundary representing the spring mesh. Left side: flat disk with α = 1. Right side: resulting elliptic shape for α > 1 (swelling of interior disk) and hyperbolic shape for α < 1 (shrinking of interior disk). Swollen springs are shown in red, while shrunk springs are shown in blue.

2.1.2 Control parameters. After all, our system of the swelling elastic disk is defined by a small set of dimensionless control parameters. These are the previously mentioned swelling factor α and the ratio of the inner and outer radius R_in/R_out. In addition, we also want to be able to describe a disk where the inner disk and the outer annulus consist of different materials.²⁵ Therefore, we introduce different elastic moduli and thus different spring constants. The spring constant k_in is valid for interior springs with r ≤ R_in, while k_out belongs to outer springs with r > R_in, and the ratio k_in/k_out is another control parameter. Finally, the thickness of the disk has an influence, even in a two-dimensional model: the relative importance of the bending energy (8) is governed by κ_B/Y_2D ∝ h², i.e., a thicker disk is harder to bend. This is usually captured by a dimensionless Föppl–von Kármán number, a dimensionless ratio of Young modulus and bending modulus, which we define for our disk as

The Föppl–von Kármán number is large for thin disks and is the fourth and last control parameter of our system.

2.2 Results

Starting with a flat disk with α = 1, we increase/decrease α in small steps Δα and minimize the energy after each step. Fig. 5(a) shows the resulting energies: the total energy, the spring energy and the bending energy (separated in mean curvature and Gaussian curvature part), as functions of α, while Fig. 5(b) shows the total height Δz of the shape. The shape can deform both into positive and negative z-direction with equal probability; we count the height Δz of the shape always as the positive absolute value of the maximal difference of z-coordinates. For small changes of α the disk stays flat at first, only the spring energy increases quadratically because of the change of the springs' rest lengths. We have E = E_s in this regime. Swelling (α > 1) or shrinking (α < 1) of the interior disk imparts elastic compression or stretching energy to the flat state, which is released in the snapping transition. At a critical swelling factor α_c2,e for increasing α (or α_c2,h for decreasing α, respectively) a transition into a curved conformation with Δz > 0 occurs. We find two stable curved configurations: for increasing α above α_c2,e > 1 the disk snaps into an elliptic (subscript “e”) dome-like shape, while it snaps into a hyperbolic (subscript “h”) saddle for decreasing α beyond α_c2,h < 1 (see Fig. 4). At these transitions, E_s is reduced, because the springs can relax to a certain degree. On the other hand, E_B is increased because of the increased curvatures in the dome- or saddle-like shapes. Increasing (decreasing) α again in order to get back to α = 1, we do not see a transition back into the flat state at α_c2,e (or α_c2,h, respectively). Instead, the shape remains curved for α < α_c2,e (α > α_c2,h). In the following, the curved disk flattens continuously with E_B and Δz decreasing until α = α_c,e (or α = α_c,h) is reached. There, E_B and Δz vanish continuously, and the disk is flat again. In conclusion, we find an apparent hysteresis loop in the deformation behavior within the red areas between α_c,e and α_c2,e (or α_c,h and α_c2,h).

Fig. 5 Energies (a), shape's height Δz and negative lowest Hessian eigenvalue −λ_min (b) in the spring mesh model as functions of the stretch factor α. Circles denote numerical values calculated with decreasing α, while crosses are related to increasing α. The disk is always flat if it is located in the white area and always curved in the cyan regions. The blue shapes illustrate the corresponding conformations of the disk. The red areas mark the regions of pseudo-hysteretic effects. Arrows illustrate the directions inside the hysteresis loops. The simulated disk had a mesh with n_B = 120 boundary vertices, R_in/R_out = 0.5, γ_FvK = 600, k_in/k_out = 0.24 and maximum fluctuations of δ_max = 5 × 10⁻⁴R_out.

2.2.1 Pseudo-hysteresis and long-wavelength bifurcation. The stability of the disk's conformation upon approaching the transition can be analyzed in more detail with the help of the eigenvalues of the Hessian matrix of the system's total energy. If the smallest eigenvalue λ_min becomes negative, there is a deformation mode leading directly to a lower energy, and the system becomes unstable. The blue scale on the right side of Fig. 5(b) shows the smallest eigenvalue of our elastic system (please note that we show the negative eigenvalue −λ_min). It is zero in the flat configuration until α exceeds α_c,e (or α_c,h). Then, still in the flat configuration in the red area, λ_min becomes significantly negative indicating that the system is unstable. This means that already in the entire red area, there is an unstable deformation mode available that leads directly into the curved conformation. After the transition, the curved conformation remains stable. Therefore, we conclude that the red area is not an area of genuine hysteresis. The transition to the curved shape could directly happen at α_c,e (or α_c,h) if the system finds the existing unstable deformation mode. The disk remains flat in the red area only for numerical reasons, and the values of α_c2,e and α_c2,h and, thus, the size of the red region actually shrinks if we increase random displacements |δ_i| (numerical fluctuations) that are imposed. In the experimental system, we expect that thermal fluctuations will always allow the disk to find the unstable mode such that hysteresis will be absent.

We can perform a linear stability analysis of the flat compressed state of the inner disk in order to further characterize the bifurcation into an elliptic dome for increasing α above α_c,e. In the limit of a small stiff outer annulus (k_in/k_out ≪ 1 and R_in/R_out ≈ 1), the effect of swelling the interior with a factor α > 1 is to establish a compressive homogeneous pre-stress σ_xx = σ_yy = −σ₀ = −Y_2D(α − 1) in the interior. We can perform a linear stability analysis of the flat state z(x,y) = 0 of an infinite plate under pre-stress −σ₀ using plate theory. Expanding the Airy stress function χ(x,y) = −σ₀(x² + y²)/2 + χ₁(x,y) and the normal displacement z(x,y) = z₁(x,y) around the flat, homogeneously pre-stressed state we find the following plate equations to linear order in χ₁ and z₁:

κ_B∇⁴z₁ + σ₀∇²z₁ = 0, ∇⁴χ₁ = 0. (15)

An Ansatz and for an oscillatory instability of the flat state with a two-dimensional wave vector q⃑ = (q_x,q_y) leads to the condition

σ₀ = Y_2D(α − 1) = κ_Bq², (16)

which is fulfilled for α > α_c,e with α_c,e − 1 ≈ κ_Bq²/Y_2D. The resulting instability is a long-wavelength instability, i.e., sets in at the smallest available wave vector q, as opposed to buckling of a spherical shell under pressure, where the pressure also leads to a homogeneous compressive pre-stress, but the buckling instability is a short-wavelength instability because of the non-vanishing background curvature.^29,38,39 For an inner disk of radius R_in ≈ R_out the shortest available wave vectors have q ∼ 1/R_out. Closer inspection shows that the unstable radially symmetric, oscillating modes are (with the Bessel function J₀). The approximate boundary condition ∂_rz₁(R_out) = 0 for a small stiff outer annulus leads towhere 3.83 is the first zero of the Bessel function J₁(x). This is in good agreement with numerical results even for k_in/k_out < 1 and R_in/R_out = 0.5 as shown in Fig. 6.

Fig. 6 Critical swelling factor α_c,e as a function of γ_FvK for R_in/R_out = 0.5 and different values of k_in/k_out. The solid black line represents the theory curve given by eqn (17).

The stability analysis with the stability eqn (15) also shows that a genuine hysteresis should be absent, and the bifurcation is a supercritical pitchfork bifurcation similar to the Euler buckling bifurcation of a beam. This is in contrast to buckling of a spherical shell under pressure, which is a subcritical bifurcation.²⁹

3 Re-establishing hysteresis

3.1 Framing the disk and additional attractive interaction

Now we want to modify the system in a way that a genuine hysteresis is re-established, which will also be present in experiments. The basic idea is to energetically penalize slightly deformed intermediate states of the disk during the transition to the curved shape resulting in an additional energy barrier for this transition. This barrier has to be overcome or decreased by additional swelling before the disk can snap into a curved state and stabilizes the flat disk.

In order to realize that, the first step is to “frame the disk”: we combine our flat disk with radius R_out in the xy-plane (Fig. 7(a)) with an additional fixed, undeformable frame (Fig. 7(b)). As a result, the boundary of the disk is now fixed (Fig. 7(c)). Therefore, the piecewise constant swelling function (4) is replaced by a globally constant, homogeneous swelling factor α for the whole disk (but not for the frame); α > 1 still corresponds to swelling the interior of the disk and, accordingly, leads to a transition into a dome-like shape (Fig. 7(d)). A saddle shape can no longer be realized in this set-up. Framing is equivalent to the above limit of a very thin and stiff outer annulus.

Fig. 7 Illustration of the set-up of a flat disk inside of a fixed frame. The flat disk (a) is placed inside a fixed and undeformable frame (b) with an attractive central region (in red). Disk and frame are compatible in the flat state (c). If the disk swells uniformly, an elliptic dome-like shape results (d).

The second step in order to create a genuine hysteresis is to introduce an additional attractive interaction between the frame and the disk. The central region of the frame A_c (red area in Fig. 7) attracts the disk leading to an additional potential energy for the disk. Inspired by an attractive van der Waals force, we choose the attractive part of a Lennard-Jones potential for the interaction,

with the total attractive potential energy or for the mesh model of the disk. The attractive potential has a finite range σ, and a potential depth −ε at z = 0; the force of this potential vanishes in the completely flat state (z = 0) improving the numerical stability. For σ ≪ R_out, on the other hand, the force nearly vanishes also in the completely deformed state because most parts of the disk are out of the potential range. In conclusion, only the transition itself is energetically penalized.

3.2 Hysteresis and short-wavelength bifurcation

In order to gain insight into the influence of an attractive potential v_pot(z) onto the instability of the swelling disk, we can consider the case where the attractive potential acts over the whole area of the disk, i.e., A_c = A. Then the linear stability analysis leads to a short-wavelength instability as eqn (15) becomes modified to

κ_B∇⁴z₁ + σ₀∇²z₁ + v_pot′′(0)z₁ = 0 (19)

resulting in an instability condition κ_Bq⁴ + σ₀q² + v_pot′′(0) < 0 (if v_pot′(0) = 0 and with v_pot′′(0) = 36 × 2^2/3ε/σ² for the potential (18)). This is exactly equivalent to the wrinkling condition of a membrane on an elastic substrate (or a Winkler foundation) under compressive stress,⁴⁰ where v_pot′′(0) corresponds to the substrate stiffness. Interestingly, this is also equivalent to the short-wavelength instability condition for buckling of a spherical shell under pressure with the homogeneous compressive pre-stress playing the role of the pre-stress from homogeneous pressure and the curvature of the potential v_pot′′(0) playing the role of the background curvature term.³⁹ Now, an instability sets in at the smallest σ₀, for which the instability condition can by fulfilled, which is for or for α > α_c,f with (subscript “f” for framed disk) and at the wave vector q₀ = (σ₀/2κ_B)^1/2 = (v_pot′′(0)/κ_B)^1/4. This is a short-wavelength instability with q₀ > 1/L if v_pot′′(0) is sufficiently large. We also expect to find a subcritical bifurcation with hysteresis in analogy to buckling of a spherical shell under pressure.²⁹

For a localized potential, i.e., if the attractive region A_c is smaller than A as in Fig. 7, we expect that the critical swelling factor is further increased such that the unstable wavelength 1/q₀ = (2κ_B/σ₀)^1/2 fits into the size of the attractive region. This results in a condition .

Analogously to Fig. 5, the behavior of the system including frame and attractive interaction is shown in Fig. 8. Increasing α starting at α = 1, the framed disk again stays flat at first until α_c2,f is reached, where the transition into the curved, dome-like shape occurs. There, we see a significant reduction of the spring energy and an increase of the bending energy. In contrast to Fig. 5, also the total energy is reduced drastically during the transition. Decreasing α again, the behavior is qualitatively the same as before, the shape continuously flattens but stays curved until α_c3,f is reached, where the disk is flat again. The significant difference to the simple set-up from above can be found by taking a look at the smallest Hessian eigenvalue λ_min (Fig. 8(b), blue scale). Between α_c3,f and α_c,f (yellow area), there are no negative eigenvalues, which means that both the flat disk and the curved shape are (meta-)stable in this region. Only if α exceeds α_c,f, λ_min becomes negative signalling an unstable flat shape. In conclusion, the region between α_c,f and α_c2,f (red area) is again a region of pseudo-hysteresis, where we see hysteresis in the numerics but hysteresis vanishes in the presence of sufficient random fluctuations and in the experiment, while the yellow area is related to a genuine hysteresis that should also be robustly observable in an experiment in the presence of some fluctuations. This gives also rise to a shape hysteresis as indicated in Fig. 1. The intermediate states upon snapping into the dome-like shape feature a flattened region around the center, while this feature is missing when the shape continuously flattens.

Fig. 8 Energies (a), shape's height Δz and negative lowest Hessian eigenvalue −λ_min (b) in the spring mesh of a disk in a fixed frame with attractive potential as functions of the swelling factor α. Crosses denote numerical values calculated with increasing α, while circles are related to decreasing α. The region with a genuine hysteresis is marked in yellow, while pseudo-hysteretic effects are again marked in red. The disk is always flat if it is located in the white area and always curved in the cyan region. Arrows illustrate the directions inside the hysteresis loops. The simulated system is the same system from Fig. 5 but with an additional potential energy for each vertex given by eqn (18). The simulated disk had a mesh with n_B = 120 boundary vertices, γ_FvK = 1066 and maximum fluctuations of δ_max = 10⁻⁴R_out. The parameters of the attractive potential were set to ε = 7.3 × 10⁻⁷k_in and σ = 0.01R_out. The potential acted in an inner region A_c with a radius of 0.2R_out.

4 Hydrodynamics

4.1 Model

In the following, we want to show that the hysteretic shape transition of the modified framed elastic disk can be exploited as a propulsion mechanism for a microswimmer under a periodic time-reversible driving of the swelling factor α. To this end, we need to model the hydrodynamic interaction between the elastic disk and a surrounding fluid. For this proof of concept we simulate the Stokesian dynamics⁴¹ and use the Rotne–Prager interaction.^42,43 This interaction describes the movement of a small sphere in the flow field of another sphere. Therefore we model our disk as a sheet of small spheres and place spheres of radius a ≪ R_out on every vertex of the spring mesh, see Fig. 9. The velocity v⃑_i of every sphere i can then be calculated from the knowledge of the external forces f⃑_j on all spheres j viaThe constant η describes the viscosity of the surrounding fluid and represents the three-dimensional unit matrix. In this model, we ignore torques and rotations for simplicity.

Fig. 9 Illustration of the sphere mesh. Small spheres (blue circles) are placed on the vertices of a Delaunay triangulation. Example mesh with n_B = 40 boundary vertices. The spheres have a radius a = 2πR_out/(10n_B), which is about 10% of their average distance.

The forces f⃑_i can be calculated (analytically) from the discretized stretching and bending energies (9) and (13) as gradients with respect to the vertex position r⃑_i at each time step, and eqn (20) gives the resulting vertex velocities. The trajectory of each sphere and the disk's center of mass as the average of all sphere positions are calculated by a simple Euler integration of the velocities with a time step Δt.

4.2 Simulation

To simulate the movement of the disk, the trajectories of the spheres are calculated based on the forces acting on them. For a fixed swelling factor α this dynamics will relax into the same force-free equilibrium state that we determined also by “dry” or static energy minimization in Fig. 8. Using the dynamics (20) we can, however, obtain a realistic dynamics of each sphere position and, thus, of the deformation and propulsion dynamics of the whole disk in the presence of hydrodynamic interactions in a viscous fluid.

The general concept of the simulation stays basically the same. Again, α is changed in small steps Δα. After each step, the trajectories of the spheres are calculated until the forces on the spheres fall below a threshold, . The simulation gives, in principle, the corresponding hydrodynamic time scale Δτ_h on which this elastic relaxation happens. That means that there are actually two different time scales operating in this system. The swelling time scale T_sw is defined by the swelling process and is the time that the disk needs to run through a complete deformation cycle. This is the time scale that can be externally controlled in an experiment, where swelling frequencies f_sw = 1/T_sw ∼ 5 Hz are possible for disks made from thermoresponsive hydrogels if plasmonic heating of embedded gold particles by laser light is utilized.³¹ If we divide a deformation cycle into N small changes Δα_n and Δτ_h,n is the hydrodynamic relaxation time for each step, we obtain the second time scale , which is the hydrodynamic relaxation time scale for one deformation cycle and determined by the interplay of elastic forces and hydrodynamic friction. In our simulation model we assume that hydrodynamic relaxation is much faster than the swelling process, τ_h ≪ T_sw, i.e., the disk swells slowly compared to its deformation motion caused by the swelling, and we can use quasi-equilibrated forces along the swelling cycle.

We can estimate an order of magnitude for the hydrodynamic relaxation time scale. The typical total force onto a disk with shape height Δz close to the instability is F ∼ Δzσ₀ (see eqn (15), which equates areal force densities); the disk has a friction coefficient ∼ηR_out, such that the typical velocity is ∂_tΔz ∼ F/ηR_out ∼ Δzσ₀/ηR_out, which leads to relaxation times

(using σ₀ ∼ Y_2D(α − 1) ∼ Y_2D/γ_FvK close to the instability for an unframed disk, see eqn (17)). Typical elastic moduli for thermoresponsive hydrogels are Y_3D ∼ 10–100 kPa⁴⁴ and Föppl–von Kármán numbers for the disks in ref. 31 are γ_FvK ≈ 300 (for R_out = 30 μm and h = 5 μm) and result in fast hydrodynamic relaxation time scales τ_h ∼ 10⁻⁵ s in water, such that swelling frequencies up to f_sw ∼ 10⁵ Hz still satisfy quasi force–equilibrium along the swelling cycle as assumed in our simulation.

On the other hand, the hydrodynamical relaxation time scale τ_h should be large enough (the hydrodynamically damped deformation or snapping velocity of the disk slow enough) for the underlying assumption of low Reynolds number hydrodynamics to apply. This is the case if Re ∼ ρR_out²/ητ_h < 1 or τ_h > ρR_out²/η. Inserting τ_h from eqn (21) we obtain a condition on the disk geometry, h³/R_out < η²/Y_3Dρ ∼ 10⁻¹ μm², where the last estimate is for the density and viscosity of water and moduli Y_3D ∼ 10 kPa. This implies that disks have to be designed sufficiently thin (and, thus, bendable) to remain at low Reynolds numbers, which is possibly a critical point for experimental realizations. For disks of radius R_out = 30 μm, thicknesses h < 2 μm are required. As long as the low Reynolds number assumption applies, the swimming distance per deformation cycle is independent of the time scale T_sw of the swelling process and, thus, the deformation velocity. The speed of shape changes affects the swimming velocity but does not affect the swimming distance as long as shape changes remain sufficiently slow that the low Reynolds number assumption holds.

The quality of our discretization and Stokesian dynamics simulation scheme can be assessed by monitoring the fluid flow

through the discretized surface (with unit normals n⃑_i), which is given by the relative velocity of the fluid flow u⃑(r⃑) with respect to the disk vertex velocities v⃑_i. The fluid flow field can be obtained from the Rotne–Prager interaction aswhile the velocity v⃑_i of vertex i is given by (20) and is essentially u⃑(r⃑_i) with a regularization of the j = i term. The quality of the approximation can be measured by the dimensionless permeability |j_S/j_P|, where j_P is the fluid stream through a theoretically perfectly permeable surface that moves with a velocity v⃑_i through the fluid (i.e., setting u⃑(r⃑) = 0 in eqn (22)). Small permeabilities indicate that discretization and Stokesian dynamics is a good approximation; ideally, we reach |j_S/j_P| = 0 because no fluid should pass through the surface. For a ≈ 0.1l₀ (where l₀ = 2πR_out/n_B is the typical rest length of a spring in the discretized mesh) we find surprisingly small permeabilities |j_S/j_P| < 20% for discretizations with n_B > 100 in view of the fact that less than 5% of the surface area are covered by spheres.

4.3 Swimming motion of the snapping elastic disk

In order to quantify the swimming motion of the disk and proof the concept of a net swimming motion, we measure the movement of the disk's center of mass as a function of time over multiple swelling cycles. Because of its symmetry, the disk only moves into the direction perpendicular to its initial plane, the z-direction. Therefore, Fig. 10 shows the z-coordinate of the center of mass, z_CoM, as a function of time for ten full swelling cycles. We will assume in the following that the dome-like shape always snaps downwards, i.e., the opening is in negative z-direction as shown in Fig. 7. In each swelling cycle α(t) the swelling factor α changes between α_min = 1 and α_max = 1.1 > α_c2,f in 200 steps back and forth in a completely time-reversible fashion.

Fig. 10 z-Coordinate z_CoM of the disk's center of mass as a function of time for ten full deformation cycles. The simulated disk featured the same parameters as the disk in Fig. 8 but with n_B = 60 and ε = 5 × 10⁻⁸k_outR_out². The spheres in the mesh had a size of a = 2πR_out/(10n_B). In each swelling cycle, the swelling factor α varied between 1 and 1.1 in 200 steps. The lower graph shows a zoom into a single (the first) cycle. The net swimming motion is in the direction of the opening.

Each swimming cycle consists of three phases (see Fig. 10). Beginning at α = 1 and z_CoM(t = 0) = 0, the disk does not move at first as long as the disk stays flat until α > α_f,c2, where the disk snaps into a dome-like shape in the numerics. The second phase starts with the short snapping process itself. During snapping into the dome-like shape, the swimmer moves into the positive z-direction (the direction of the tip of the elliptic dome, upwards in Fig. 7). This movement happens nearly instantaneously on the swelling time scale. Then, for α_f,c2 < α < α_max, the height Δz of the dome-shape slowly continues to grow. Then, the swelling factor starts to shrink again, and the third phase starts. Here, the shape flattens again, and we see a movement into the negative z-direction. When the disk arrives back in its original state (α = 1), a small net motion into the negative z-direction, i.e. in the direction of the opening, remains, which results in an effective slope in the diagram over multiple swelling cycles in Fig. 10.

The effective speed of this swimmer is quite small with a net swimming distance of just below 1% of the disks radius R_out after ten deformation cycles for the parameter values given in Fig. 10. This is a consequence of the known problem of all shape-changing microswimmers that the swimming distance typically scales only quadratically with the deformation displacement.^11,12 Nevertheless, this proves the principle that a the hysteretic shape transition of a flat elastic disk can be used as a propulsion mechanism for a microswimmer. Shape-changing swimmers can still be effective if the driving frequency (the swelling cycle frequency for our swimmer) is sufficiently high. At low Reynolds numbers, the resulting swimming distance Δz_CoM = z_CoM(t = T_sw) is also independent of the shape of the actual swelling cycle α(t) but can only depend on the values of α_min = 1 and α_max = α_f,c2 characterizing the hysteretic part of the swelling cycle.

We can use the reciprocal theorem⁴⁵ in order to obtain analytical insight into the dependence of the swimming distance Δz_CoM on the deformation displacement Δz of the swimmer. We focus on the z-component as axisymmetric deformations can only lead to motion along the axis of symmetry of the disk. For a deformed disk shape A(t) with a shape deformation z = z(r,t) and a deformation velocity ż in z-direction (r is the radial coordinate in the xy-plane) the reciprocal theorem gives

where F_D is the z-component of the viscous drag force of a disk with rigid shape A(t) and velocity v_CoM in z-direction, σ_D the corresponding stress tensor of the fluid, and (n⃑·σ)_z the z-component of the normal stress of the fluid onto the shape A(t); these quantities are related via. This leads towhereis the stress difference (divided by drag force) between the two faces of a disk with rigid shape A(t). The first term in the last equation is the result for the flat disk (see ref. 46–48). For weakly deformed rigid disks there are corrections; for the scaling of the net swimming distance with the shape height Δz it is crucial how these corrections scale with Δz. For the total drag force the corrections are quadratic, −O(F_D(z/R_out)²). This can be shown by interpreting the weakly deformed rigid disk as perturbed flat disk and applying the reciprocal theorem.⁴⁹ Because both the fluid velocity field and the pressure only vary quadratically in z close to the flat disk shape,⁴⁷ boundary conditions to the flow and, thus, the fluid stresses and the drag receive only quadratic corrections. This is also supported by an exact result for the axisymmetric Stokes flow past spherical caps of opening angle β and radius R,⁵⁰ which can also be interpreted as weakly bent rigid disks. The friction force F_D,cap = μv_CoMR(6β + 8 sin β + sin(2β)) reduces in the disk limit β ≪ 1, where R_out ≈ Rβ and z/R_out ≈ β/2, to F_D ≈ 16μv_CoMR_out(1 − (z/R_out)²/6) in the first two leading orders. Therefore, we expect to find a leading order reduction of the stress jump (26), which is quadratic in z/R_out for a weakly deformed rigid disk. This allows us to extract the scaling of the swimming velocity as a function of the deformation function z(r,t) which describes the shape hysteresis.

The shape hysteresis as sketched in Fig. 1 is mainly caused by flattened shapes during snapping (ż(r) < 0, v_CoM > 0), while the shape remains dome-like during re-flattening (ż(r) > 0, v_CoM < 0). Expanding in the deformation z in eqn (25) and integrating over time to obtain the swimming distance, , we realize that the leading order term only depends on final and initial state. It gives the same swimming distance Δz_CoM,snap ∼ −Δz both in the snapping and re-flattening phase and cancels out exactly for a full shape cycle. The next-to-leading term, however, gives a slightly bigger contribution in the re-flattening phase, which should scale as Δz_CoM/R_out ∼ −(Δz/R_out)³. Therefore, the net swimming distance is only a third order effect in the hysteretic shape height Δz. This is even smaller than the quadratic order observed, for example, for deforming spheres^11,12 and appears to be a hydrodynamic consequence of the disk geometry. This is confirmed in Fig. 11, where the net swimming distance for swimmers with different thicknesses is shown as a function of their shape height in the snapped state.

Fig. 11 Net swimming distance |Δz_CoM/R_out| per cycle as a function of the shape height Δz/R_out both for disks and 9-bead swimmer. For the 9-bead swimmer the parameter Δz is directly prescribed. For disks the increasing shape height Δz in the snapped state is realized by increasing the thickness in the range h/R_out = 0.02 − 0.14, the other parameters are as in Fig. 8. The fit is Δz_CoM/R_out = 0.022(Δz/R_out)³.

The shape height in the subcritical bifurcation with hysteresis is limited by the stretching factor α of the shape. Because a fiber along the diameter of the disk will stretch by a factor α − 1, the change in height will scale as Δz/R_out ∼ (α − 1)^1/2 for a curved disk. For the hysteretic part of the deformation cycle, the relevant stretching factor is α = α_c,f. This results in a net swimming distance

The swimming distance Δz_CoM per stroke can be increased by increasing α_c,f − 1 and, thus, the width of the yellow hysteresis area in Fig. 8 which can be achieved, for example, by increasing the thickness and, thus, the bending modulus of the disk. This is explicitly demonstrated in Fig. 11. The thickness dependence α_c,f − 1 ∝ h² results in |Δz_CoM| ∝ h³.

The deformation of the disk into a dome-like shape resembles the deformation pattern of a scallop as envisioned by Purcell;² our disk is, however, hysteretic, which enables swimming. Interestingly, we find swimming into negative z-direction, which is the direction away from the tip of the dome. This is the opposite swimming direction that we expect from Purcells's high Reynolds number scallop, which moves in the direction of its tip by thrusting fluid out of the opening. Our disk, on the contrary, will pull in fluid through the opening because the size of the opening is fixed by the frame. Also at high Reynolds numbers this leads to inertial thrust in the direction of the opening.

5 9-bead model

In order to further elucidate the underlying propulsion mechanism, we present a strongly simplified 9-bead model of the hysteretic scallop. In this model, the swimmer only consists of nine spheres that move on simple prescribed trajectories in relative coordinates. The 9-bead-model features three deformation phases, which mimic the above three phases of the disk swimmer; these deformation phases are visualized in Fig. 12. The sketches in the first row show the corresponding states of the full, framed disk analogously to Fig. 7. The second row shows three-dimensional illustrations of the 9-bead swimmer and the last row a two-dimensional side view.

Fig. 12 Illustration of the cyclic deformation in the 9-bead model. First row: corresponding states of the disk in the fixed frame. Second row: three-dimensional sketch of the 9-bead swimmer sates. Third row: two-dimensional side view on the 9-bead model. Each column represents one deformation state (1–2–3–1).

The undeformed flat state is represented in the left column. Nine spheres are aligned symmetrically in the xy-plane. The central sphere (red) is surrounded by four spheres (blue) that are placed in a symmetrical way in a distance of R_out/2. The last four spheres (green) are set in an identical way but in a distance of R_out to the central sphere. They represent the fixed frame. The connection lines between the spheres indicate the structure of the object without a direct physical meaning. The second state (second column) represents a typical transition state of the disk during the snapping transition into the dome-like shape. There, the outer regions of the disk are already close to their final position, but the central region does not show a tip yet. This is modelled by elevating the five inner spheres (red and blue) in z-direction with Δz. The third state then is the finally stable elliptic dome-like shape after the transition. Here, the central sphere is again elevated by another distance of Δz. Finally, the last state is again the flat and relaxed disk so that the cycle can be repeated. During the evolution from one state into another, the spheres move on trajectories, which simply linearly interpolate the positions of the spheres between initial and final state. Therefore, the first evolution phase (1 → 2) and the second phase (2 → 3) represent an increasing swelling factor α, while the last phase (3 → 1) is related to a decreasing α.

5.1 Swimming behavior

The swimming motion of the 9-bead model is calculated with the same hydrodynamic model as before. But now, the relative velocities of the spheres are prescribed. So the forces on them have to be calculated. The knowledge of eight relative velocities, together with the condition that the swimmer is force-free, i.e., the total force must vanish, enables us to calculate the motion of the center of mass.^17,18 The results are shown in Fig. 13 for a single deformation cycle. In the first deformation phase (0 < t < 1/3T_sw), the swimmer moves upwards in positive z-direction. In the second phase, when the central sphere forms the tip (1/3T_sw < t < 2/3T_sw), there is still an upwards movement, which is simply slower than before because less spheres show a relative motion in this phase. In the final phase (t > 2/3T_sw), all spheres relax back to their original position causing a quick movement in negative direction. After one complete cycle, at t = T_sw, we see a small net movement in negative direction with Δz_CoM = −6.6 × 10⁻⁴R_out. In conclusion, we see the same qualitative behavior as for the full disk model in Fig. 10. We can state that the simplified 9-bead model is able to model the swimming mechanism of a swelling disk that has a fixed outer frame. For the 9-bead swimmer simulation results in Fig. 11 show that the swimming distance per cycle is also compatible with z_CoM/R_out ∝ −(Δz/R_out)³ for small prescribed shape heights Δz and collapses onto the results for the swelling disk.

Fig. 13 z-Coordinate of the 9-bead swimmer's center of mass as a function of time during a single deformation cycle. The bead size is a = 0.1R_out and Δz = 0.5R_out. The sketches illustrate the conformations of the swimmer during the deformation cycle.

We can also compare the characteristics of the resulting fluid flow field between the full disk model and the 9-bead model; Fig. 14 shows that there is good qualitative agreement between both models. In particular, the characteristic stagnation point and the vortex ring below the snapping disk are reproduced. These features are characteristic for non-convex bodies and are also observed for dragged spherical caps.⁵⁰

Fig. 14 Velocity field u⃑(r⃑) in the (y = 0) -plane for each deformation phase. The left side shows the disk with a fixed frame from Fig. 8 and a discretization of n_B = 40. The right side shows a 9-bead swimmer with a = 0.1R_out and Δz = 0.5R_out. Colors indicate the absolute |u⃑| in units of k_out/η for the disk swimmer and in units of R_out/T_sw for the 9-bead swimmer. The arrows show the direction of the fluid velocity vectors.

We can expand the axisymmetric flow field in the far-field into Legendre polynomials^11,12 and extract the dipole contribution p₂P₂(cos θ)/r² of the radial part u_r(r,θ) of the flow field u⃑(r⃑) from eqn (23) and normalize by the center of mass velocity |v_CoM| to define the dimensionless parameter

Values γ < 0 (p₂ < 0) indicate “pulling” motion of the swimmer, values γ > 0 (p₂ < 0) “pushing” motion, and γ ≈ 0 (p₂ ≈ 0) a neutral motion.³ The results for γ for both 9-bead and disk swimmer are shown in Table 1. In the beginning of phase 1 → 2 (state 1 in Fig. 12) the disk is approximately neutral, while the 9-bead swimmer is a weak pusher; the disk appears neutral because the subcritical snapping instability is caused by short-wavelength deformations. In the beginning of phase 2 → 3 (state 2 in Fig. 12) disk and 9-bead swimmer are pullers, and in the beginning of phase 3 → 1 (state 3 in Fig. 12) they are pushers.

Table 1 Swimmer parameter γ (see eqn (28)) for 9-bead and disk swimmer in the beginning of the respective phases

Phase	γ 9-bead	γ disk
1 → 2 (neutral/weakly pushing, γ small)	−0.0068	0.39
2 → 3 (pulling, γ < 0)	−3.63	−10.88
3 → 1 (pushing, γ > 0)	9.18	1.24

---

6 Conclusions

We pursued the concept of a low Reynolds number swimmer based on the concept of utilizing a hysteretic shape transition in order to convert a completely time-reversible oscillation of a control parameter into a directed swimming motion. We proved this concept with a flat circular elastic disk that undergoes a shape transition into curved shapes by a localized swelling of an inner disk or an exterior annulus. While swelling the inner region of the disk with a constant swelling factor α and the outer annulus with 1/α, we saw a transition into an elliptic dome-like shape for α > 1 and a transition into a hyperbolic saddle shape for α < 1. The control parameter of this shape transition is the swelling factor α. We found the transition to be a supercritical bifurcation with only numerical hysteresis, which will disappear in a real experiment in the presence of some fluctuations. We could re-establish a genuinely hysteretic shape transition by replacing the outer annulus by a fixed outer frame and by introducing an additional attractive short-range interaction in the central region. The details of this interaction are not relevant, as long as it is an attractive short-range interaction, as it arises, for example, from van der Waals forces or screened electrostatic forces. We then see a hysteretic subcritical shape transition between a flat state and a dome-like shape.

Embedding this framed disk into a viscous fluid at low Reynolds numbers, a Stokesian dynamics simulation of the hydrodynamic interaction with the surrounding fluid with a time-reversible and cyclic changing of α showed that the swimmer is effectively moving into the direction of the opening of the dome. This way, a self-deforming microswimmer can be realized that uses only a single scalar control parameter, the swelling factor α. This control parameter has to be changed only by a few percent in order to trigger a drastic conformation change with the snapping transition into the curved shape. Interestingly, the snapping into an elliptic dome-like shape resembles the opening and closing of a scallop as envisioned in the scallop theorem.² As opposed to Purcell's scallop the elastic disk swimmer actually performs directed motion at low Reynolds numbers because of the additional hysteresis. The swimming direction of the snapping elastic disk is into the direction of the opening of the dome (away from the tip). As for many shape-changing low Reynolds number microswimmers swimming is a higher order effect in the deformation displacement. For the snapping elastic disk the net swimming distance per swelling cycle is only a third order effect in the height of the dome-like shape (|Δz_CoM/R_out| ∼ −(Δz/R_out)³). The swimming mechanism by hysteretic snapping is reproduced by a simplified 9-bead model of the disk. Qualitative agreement of the resulting flow fields, its pusher/puller characteristics and the net swimming distance (see Fig. 11 and 14) show that the 9-bead model captures the essence of the swimming mechanism.

An experimental realization of the snapping disk microswimmer concept could be possible using thermoresponsive hydrogels, which have already been used for the implementation of helical microswimmers.^31,32 These hydrogels can be swollen by plasmonic heating of embedded gold particles by laser light. Therefore, localized swelling is possible by embedding gold particles only in specific parts of the hydrogel, such that these hydrogels are suited to realize the proposed disk microswimmers. With hydrogel properties from the literature, an estimate of the expected swimming speed is also possible. Mourran et al. investigated the thermal swelling behavior of thermoresponsive hydrogels with the help of circular disks.³¹ These disks had a diameter of 30 μm and a thickness of 5 μm. They showed cyclic diameter changes of several percent caused by a pulsed laser with frequencies of 1–5 Hz. Therefore, a frequency of 5 Hz seems to be realistic for the realization of the deformation cycle of our microswimmer with comparable dimensions. With a diameter of 30 μm and a swimming distance of about 1% of the radius after ten cycles, a net velocity of 0.075 μm s⁻¹ follows. Generally, shape-changing swimmers such as the snapping disk can only be made effective if the driving frequency (the swelling cycle frequency for our swimmer) is sufficiently high. On the one hand, we expect higher frequencies to be possible for thin disks. On the other hand, the swimming distance per swelling cycle can be increased by increasing the thickness of the disk (|Δz_CoM| ∝ h³).

One has to bear in mind that the snapping transition dynamics has to be sufficiently slow for low Reynolds number hydrodynamics to apply. We estimated the typical hydrodynamical relaxation time scale, which is the typical time scale for snapping in a viscous fluid, to be τ_h ∼ ηR_out/σ₀ ∼ γ_FvK^3/2η/Y_3D, which is of the order of τ_h ∼ 10⁻⁵ s for the thermoresponsive hydrogel disks from ref. 31. Then Re < 1 is equivalent to τ_h > ρR_out²/η and restricts the low Reynolds number regime to sufficiently thin disks with h³/R_out < 10⁻¹ μm². At higher Reynolds numbers, we also expect swimming motion in the direction away from the tip of the dome.

The concept of a low Reynolds number swimmer based on the principle of a hysteretic shape transition can be applied to other snapping systems such as tube-like⁵¹ or shell-based snapping systems.⁵² All systems with a hysteretic snapping instability should give rise to a net swimming motion in a viscous fluid at low Reynolds numbers. As opposed to our disk swimmer, these tube and shell objects enclose volume and employ volume or pressure control, such that genuinely autonomous swimming is not possible because swimming requires attached devices to control volume or pressure.²⁶ Swelling or shrinking materials as proposed in this work are an alternative to control the equilibrium area and, consequently, also equilibrium enclosed volume of these objects by swelling and shrinking, which can often trigger the same snapping transition at fixed actual volume.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge financial support by the Deutsche Forschungsgemeinschaft via SPP 1726 “Microswimmers” (KI 662/7-2).

Notes and references

1. P. Illien, R. Golestanian and A. Sen, Chem. Soc. Rev., 2017, 46, 5508–5518.
2. E. M. Purcell, Am. J. Phys., 1977, 45, 3–11.
3. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 96601.
4. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
5. G. Taylor, Proc. R. Soc. London, Ser. A, 1951, 209, 447–461.
6. H. C. Berg and R. A. Anderson, Nature, 1973, 245, 380–382.
7. R. E. Goldstein, Annu. Rev. Fluid Mech., 2015, 47, 343–375.
8. R. Jeanneret, M. Contino and M. Polin, Eur. Phys. J.: Spec. Top., 2016, 225, 2141–2156.
9. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 862–865.
10. S. Tottori and B. J. Nelson, Biomicrofluidics, 2013, 7, 061101.
11. M. Lighthill, Commun. Pure Appl. Math., 1952, 5, 109–118.
12. J. Blake, J. Fluid Mech., 1971, 46, 199–208.
13. A. Shapere and F. Wilczek, Phys. Rev. Lett., 1987, 58, 2051–2054.
14. B. U. Felderhof and R. B. Jones, Phys. A, 1994, 202, 119–144.
15. J. E. Avron, O. Gat and O. Kenneth, Phys. Rev. Lett., 2004, 93, 1–4.
16. A. A. Evans, S. E. Spagnolie and E. Lauga, Soft Matter, 2010, 6, 1737–1747.
17. A. Najafi and R. Golestanian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 062901.
18. R. Golestanian and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 036308.
19. R. Ledesma-Aguilar, H. Löwen and J. M. Yeomans, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 70.
20. M. S. Rizvi, A. Farutin and C. Misbah, Phys. Rev. E, 2018, 97, 023102.
21. M. Leoni, J. Kotar, B. Bassetti, P. Cicuta and M. C. Lagomarsino, Soft Matter, 2009, 5, 472–476.
22. R. Golestanian, Phys. Rev. Lett., 2010, 105, 018103.
23. Y. Klein, E. Efrati and E. Sharon, Science, 2007, 315, 1116–1120.
24. E. Efrati, E. Sharon and R. Kupferman, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 016602.
25. M. Pezzulla, S. A. Shillig, P. Nardinocchi and D. P. Holmes, Soft Matter, 2015, 11, 5812–5820.
26. A. Djellouli, P. Marmottant, H. Djeridi, C. Quilliet and G. Coupier, Phys. Rev. Lett., 2017, 119, 224501.
27. S. Knoche and J. Kierfeld, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 046608.
28. S. Knoche and J. Kierfeld, Soft Matter, 2014, 10, 8358–8369.
29. L. Baumgarten and J. Kierfeld, Phys. Rev. E, 2019, 99, 022803.
30. S. Palagi, A. G. Mark, S. Y. Reigh, K. Melde, T. Qiu, H. Zeng, C. Parmeggiani, D. Martella, A. Sanchez-Castillo, N. Kapernaum, F. Giesselmann, D. S. Wiersma, E. Lauga and P. Fischer, Nat. Mater., 2016, 15, 647–653.
31. A. Mourran, H. Zhang, R. Vinokur and M. Möller, Adv. Mater., 2017, 29, 1604825.
32. H. Zhang, A. Mourran and M. Möller, Nano Lett., 2017, 17, 2010–2014.
33. L. Koens, H. Zhang, M. Moeller, A. Mourran and E. Lauga, Eur. Phys. J. E: Soft Matter Biol. Phys., 2018, 41, 119.
34. B. Kornberger, Fade2D Delauney Trianglulation library, 2016, https://protect-eu.mimecast.com/s/iX5TC0YrLfwVDBcwyCGS?domain=geom.at.
35. M. Ostoja-Starzewski, Appl. Mech. Rev., 2002, 55, 35–60.
36. H. S. Seung and D. R. Nelson, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 1005–1018.
37. K. A. Brakke, Exp. Math., 1992, 1, 141–165.
38. J. W. Hutchinson, J. Appl. Mech., 1967, 34, 49.
39. L. Baumgarten and J. Kierfeld, Phys. Rev. E, 2018, 97, 052801.
40. Z. Y. Huang, W. Hong and Z. Suo, J. Mech. Phys. Solids, 2005, 53, 2101–2118.
41. L. Durlofsky, J. F. Brady and G. Bossis, J. Fluid Mech., 1987, 180, 21–49.
42. J. Dhont, An Introduction to Dynamics of Colloids, Elsevier Science, 1996.
43. J. Rotne and S. Prager, J. Chem. Phys., 1969, 50, 4831–4837.
44. T. R. Matzelle, G. Geuskens and N. Kruse, Macromolecules, 2003, 36, 2926–2931.
45. H. Stone and A. Samuel, Phys. Rev. Lett., 1996, 77, 4102–4104.
46. S. C. Gupta, Z. Angew. Math. Phys., 1957, 8, 257–261.
47. J. P. Tanzosh and H. A. Stone, Chem. Eng. Commun., 1996, 148–150, 333–346.
48. B. U. Felderhof, J. Chem. Phys., 2012, 137, 084906.
49. H. Masoud and H. A. Stone, J. Fluid Mech., 2019, 879, P1.
50. J. M. Dorrepaal, M. E. O'neill and K. B. Ranger, J. Fluid Mech., 1976, 75, 273–286.
51. J. T. Overveldea, T. Kloeka, J. J. D'Haena and K. Bertoldia, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 10863–10868.
52. B. Gorissen, D. Melancon, N. Vasios, M. Torbati and K. Bertoldi, Sci. Robot., 2020, 5, 1967.

---