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

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.


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. 1 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 timereversal: 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. 2 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][6][7][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 2 and including swimmers performing small deformations of spheres, circles, or cylinders, [11][12][13][14][15] or shape-changing vesicles. 16 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) 19 and steer. 20 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 threebead 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][24][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][28][29] which will turn out to be conceptually similar to the elastic instability triggered in the elastic disk by localized 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. 31 These ideas were used to design swimmers with a helical structure that can be propelled by conformation changes of the helix. [31][32][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 timereversibility 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

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 twodimensional disk in polar coordinates (r,j). 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][24][25] By swelling we mean a local isotropic swelling, where the rest lengths of fibers change by a position-dependent factor A(r,j) independent of fiber orientation. In the following, 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).
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) 4 1 corresponds to local swelling, A(r) o 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 r Ð r 0 AðrÞdr, 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 j remains unchanged. In Gaussian coordinates a deformed fiber in radial direction has length dl r = dr = A(r)dr, a deformed circumferential fiber a length dl j = rA(r)dj. In general, the fiber length is related to the metric by dl 2 = g rr dr 2 + 2g rj drdj + g jj dj 2 such that the metric tensor of the swollen disk can be read off as in Gaussian normal coordinates. This is the so-called target metric which represents the preferred equilibrium state of the swollen disk. 24 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 r and j of the metric) we find rðrÞAðrðrÞÞ ; (2) as a function of the Gaussian radial coordinate r. We can transform to the standard radial coordinate r by using dr/dr = 1/A(r), which gives KðrÞ ¼ À A 0 ðrÞ þ rA 00 ðrÞ À rA 02 ðrÞ=AðrÞ rA 3 ðrÞ : 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., 25 we divide the disk into two parts: An inner disk with 0 r r r R in and an outer annulus with R in o r 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 a, while the outer annulus shall always do the exact opposite, i.e., it shrinks with the inverse constant factor 1/a. In total, the considered swelling functions A(r) can be written as ; r 2 ðR in ; R out : The swelling process is thus defined by two simple control parameters: the swelling factor a 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 a. We can distinguish between two general cases: (a) a 4 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) a o 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 a 4 1, we have A 0 (r) = 0 except around r = R in , where ÀA 0 (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 a 4 1 and negative for a o 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 a 4 1, and a hyperbolic target shape with a negative Gaussian curvature for a o 1.
Obviously, for a a 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 threedimensional embedding space and is an example of non-Euclidean geometry. 24 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 % g 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. 24 The elastic energy contains a stretching and a bending contribution. The stretching contribution is caused by strains e ij = (g ij À % g ij )/2, where g is the actual metric that the deformed swollen state assumes, and given by 24 A ijkl ¼ Y 2D 1 À n 2 n g ij g kl þ 1 À n 2 g ik g jl þ g il g jk À Á ; where we use Einstein summation, raising of indices is performed with the target metric, and (u 1 ,u 2 ) = (r,j) in Gaussian normal parametrization. The elastic tensor A ijkl is given by the two dimensional Young modulus Y 2D and the Poisson ratio n, 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 24 where the last line applies to % L = 0, H is the mean curvature, K the Gaussian curvature, and k B = Y 2D h 2 /(12(1 À n 2 )) 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 % g ij E g ij in A ijkl . Typical strains e ij are p (1 À a) 2 such that corrections are O((1 À a) 2 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 34 ), 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 twodimensional Young modulus Y 2D is given by the spring constant k and the Poisson ratio n 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: 37 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 r i A 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 y 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 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 displacementr i !r i þd i with jd i j ( l. 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.

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 a 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. 25 Therefore, we introduce different elastic moduli and thus different spring constants. The spring constant k in is valid for interior springs with r r R in , while k out belongs to outer springs with r 4 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)  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 y j and is used in the calculation of the Gaussian curvature.
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.

Results
Starting with a flat disk with a = 1, we increase/decrease a in small steps Da 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 a, while Fig. 5(b) shows the total height Dz of the shape. The shape can deform both into positive and negative z-direction with equal probability; we count the height Dz of the shape always as the positive absolute value of the maximal difference of z-coordinates. For small changes of a 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 (a 4 1) or shrinking (a o 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 a c2,e for increasing a (or a c2,h for decreasing a, respectively) a transition into a curved conformation with Dz 4 0 occurs. We find two stable curved configurations: for increasing a above a c2,e 4 1 the disk snaps into an elliptic (subscript ''e'') dome-like shape, while it snaps into a hyperbolic (subscript ''h'') saddle for decreasing a beyond a c2,h o 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) a again in order to get back to a = 1, we do not see a transition back into the flat state at a c2,e (or a c2,h , respectively). Instead, the shape remains curved for a o a c2,e (a 4 a c2,h ). In the following, the curved disk flattens continuously with E B and Dz decreasing until a = a c,e (or a = a c,h ) is reached. There, E B and Dz 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 a c,e and a c2,e (or a c,h and a c2,h ).
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 l 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 Àl min ). It is zero in the flat configuration until a exceeds a c,e (or a c,h ). Then, still in the flat configuration in the red area, l 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 a c,e (or a 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 a c2,e and a c2,h and, thus, the size of the red region actually shrinks if we increase random displacements |d 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 a above a c,e . In the limit of a small stiff outer annulus (k in /k out { 1 and R in /R out E 1), the effect of swelling the interior with a factor a 4 1 is to establish a compressive homogeneous pre-stress s xx = s yy = Às 0 = ÀY 2D (a À 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 Às 0 using plate theory. Expanding the Airy stress function w(x,y) = Às 0 (x 2 + y 2 )/2 + w 1 (x,y) and the normal displacement z(x,y) = z 1 (x,y) around the flat, homogeneously pre-stressed state we find the following plate equations to linear order in w 1 and z 1 : An Ansatz w 1 ¼ ae iqÁr and z 1 ¼ be iqÁr for an oscillatory instability of the flat state with a two-dimensional wave vector q = (q x ,q y ) leads to the condition which is fulfilled for a 4 a c,e with a c,e À 1 E k B q 2 /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 E R out the shortest available wave vectors have q B 1/R out . Closer inspection shows that the unstable radially symmetric, oscillating modes are z 1 ðrÞ ¼ bJ 0 r ffiffiffiffiffiffiffiffiffiffiffiffi ffi s 0 =k B p (with the Bessel function J 0 ). The approximate boundary condition q r z 1 (R out ) = 0 for a small stiff outer annulus leads to a c;e À 1 % 3: where 3.83 is the first zero of the Bessel function J 1 (x). This is in good agreement with numerical results even for k in /k out o 1 and R in /R out = 0.5 as shown in Fig. 6. 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. 29 3 Re-establishing hysteresis

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 a for the whole disk (but not for the frame); a 4 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.  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 E pot ¼ Ð Ac dAv pot ðzÞ or E pot ¼ P i2Ac A i v pot ðz i Þ for the mesh model of the disk. The attractive potential has a finite range s, and a potential depth Àe at z = 0; the force of this potential vanishes in the completely flat state (z = 0) improving the numerical stability. For s { 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.

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 k B r 4 z 1 + s 0 r 2 z 1 + v pot 00 (0)z 1 = 0 (19) resulting in an instability condition k B q 4 + s 0 q 2 + v pot 00 (0) o 0 (if v pot 0 (0) = 0 and with v pot 00 (0) = 36 Â 2 2/3 e/s 2 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, 40 where v pot 00 (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 00 (0) playing the role of the background curvature term. 39 Now, an instability sets in at the smallest s 0 , for which the instability condition can by fulfilled, which is for s 0 4 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k B v pot 0 0 ð0Þ p or for a 4 a c,f with This is a short-wavelength instability with q 0 4 1/L if v pot 00 (0) is sufficiently large. We also expect to find a subcritical bifurcation with hysteresis in analogy to buckling of a spherical shell under pressure. 29 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 0 = (2k B /s 0 ) 1/2 fits into the size ffiffiffiffiffi ffi A c p of the attractive region.
This results in a condition s 0 4 max 2k B .
Analogously to Fig. 5, the behavior of the system including frame and attractive interaction is shown in Fig. 8. Increasing a starting at a = 1, the framed disk again stays flat at first until a 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 a again, the behavior is qualitatively the same as before, the shape continuously flattens but stays curved until a 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 l min (Fig. 8(b), blue scale). Between a c3,f and a 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 a exceeds a c,f , l min becomes negative signalling an unstable flat shape. In conclusion, the region between a c,f and a c2,f (red area) is again a region of pseudo-hysteresis, where Fig. 8 Energies (a), shape's height Dz and negative lowest Hessian eigenvalue Àl min (b) in the spring mesh of a disk in a fixed frame with attractive potential as functions of the swelling factor a. Crosses denote numerical values calculated with increasing a, while circles are related to decreasing a. 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, g FvK = 1066 and maximum fluctuations of d max = 10 À4 R out . The parameters of the attractive potential were set to e = 7.3 Â 10 À7 k in and s = 0.01R out . The potential acted in an inner region A c with a radius of 0.2R out . 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.

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 a. 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 41 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 ðr i Àr j Þ ðr i Àr j Þ jr i Àr j j 2 f j : The constant Z describes the viscosity of the surrounding fluid and I represents the three-dimensional unit matrix. In this model, we ignore torques and rotations for simplicity.
The forces f i can be calculated (analytically) from the discretized stretching and bending energies (9) and (13)

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 a 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, a is changed in small steps Da. After each step, the trajectories of the spheres are calculated until the forces on the spheres fall below a threshold, P jf i j o e. The simulation gives, in principle, the corresponding hydrodynamic time scale Dt 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 B 5 Hz are possible for disks made from thermoresponsive hydrogels if plasmonic heating of embedded gold particles by laser light is utilized. 31 If we divide a deformation cycle into N small changes Da n and Dt h,n is the hydrodynamic relaxation time for each step, we obtain the second time scale t h ¼ P n Dt h;n , 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, t 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 P jf i j o e 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 Dz close to the instability is F B Dzs 0 (see eqn (15), which equates areal force densities); the disk has a friction coefficient BZR out , such that the typical velocity is q t Dz B F/ZR out B Dzs 0 /ZR out , which leads to relaxation times (using s 0 B Y 2D (a À 1) B Y 2D /g FvK close to the instability for an unframed disk, see eqn (17)). Typical elastic moduli for thermoresponsive hydrogels are Y 3D B 10-100 kPa 44 and Föppl-von Kármán numbers for the disks in ref. 31 are g FvK E 300 (for R out = 30 mm and h = 5 mm) and result in fast hydrodynamic relaxation time scales t h B 10 À5 s in water, such that swelling frequencies up to f sw B 10 5 Hz still satisfy quasi force-equilibrium along the swelling cycle as assumed in our simulation. On the other hand, the hydrodynamical relaxation time scale t 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 B rR out 2 /Zt h o 1 or t h 4 rR out 2 /Z. Inserting t h from eqn (21) we obtain a condition on the disk geometry, h 3 /R out o Z 2 /Y 3D r B 10 À1 mm 2 , where the last estimate is for the density and viscosity of water and moduli Y 3D B 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 mm, thicknesses h o 2 mm 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 4jr Àr j j I þ ðr Àr j Þ ðr Àr j Þ jr i Àr j j 2 þ a 3 4jr Àr j j 3 I À 3 ðr Àr j Þ ðr Àr j Þ jr Àr j j 2 f j ; while 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  (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 E 0.1l 0 (where l 0 = 2pR out /n B is the typical rest length of a spring in the discretized mesh) we find surprisingly small permeabilities | j S /j P | o 20% for discretizations with n B 4 100 in view of the fact that less than 5% of the surface area are covered by spheres.

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 a(t) the swelling factor a changes between a min = 1 and a max = 1.1 4 a c2,f in 200 steps back and forth in a completely time-reversible fashion.
Each swimming cycle consists of three phases (see Fig. 10). Beginning at a = 1 and z CoM (t = 0) = 0, the disk does not move at first as long as the disk stays flat until a 4 a 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 a f,c2 o a o a max , the height Dz 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 (a = 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 The spheres in the mesh had a size of a = 2pR out /(10n B ). In each swelling cycle, the swelling factor a 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. 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. Shapechanging 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 Dz CoM = z CoM (t = T sw ) is also independent of the shape of the actual swelling cycle a(t) but can only depend on the values of a min = 1 and a max = a f,c2 characterizing the hysteretic part of the swelling cycle.
We can use the reciprocal theorem 45 in order to obtain analytical insight into the dependence of the swimming distance Dz CoM on the deformation displacement Dz 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 : z 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, s D the corresponding stress tensor of the fluid, and ( -nÁs) z the zcomponent of the normal stress of the fluid onto the shape A(t); these quantities are related via where ½sðr; tÞ ðñ Á sÞ z F D is 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][47][48]. For weakly deformed rigid disks there are corrections; for the scaling of the net swimming distance with the shape height Dz it is crucial how these corrections scale with Dz. For the total drag force F D ¼ Ð AðtÞ dAðñ Á s D Þ z the corrections are quadratic, ÀO(F D (z/R out ) 2 ). This can be shown by interpreting the weakly deformed rigid disk as perturbed flat disk and applying the reciprocal theorem. 49 Because both the fluid velocity field and the pressure only vary quadratically in z close to the flat disk shape, 47 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 b and radius R, 50 which can also be interpreted as weakly bent rigid disks. The friction force F D,cap = mv CoM R(6b + 8 sin b + sin(2b)) reduces in the disk limit b { 1, where R out E Rb and z/R out E b/2, to F D E 16mv CoM R out (1 À (z/R out ) 2 /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 ( : z(r) o 0, v CoM 4 0), while the shape remains dome-like during re-flattening ( : z(r) 4 0, v CoM o 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 Dz CoM,snap B ÀDz 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 Dz CoM /R out B À(Dz/R out ) 3 . Therefore, the net swimming distance is only a third order effect in the hysteretic shape height Dz. 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.
The shape height in the subcritical bifurcation with hysteresis is limited by the stretching factor a of the shape. Because a fiber along the diameter of the disk will stretch by a factor a À 1, the change in height will scale as Dz/R out B (a À 1) 1/2 for a curved disk. For the hysteretic part of the deformation cycle, Fig. 11 Net swimming distance |Dz CoM /R out | per cycle as a function of the shape height Dz/R out both for disks and 9-bead swimmer. For the 9-bead swimmer the parameter Dz is directly prescribed. For disks the increasing shape height Dz 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 Dz CoM /R out = 0.022(Dz/R out ) 3 . the relevant stretching factor is a = a c,f . This results in a net swimming distance The swimming distance Dz CoM per stroke can be increased by increasing a 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 a c,f À 1 p h 2 results in |Dz CoM | p h 3 .
The deformation of the disk into a dome-like shape resembles the deformation pattern of a scallop as envisioned by Purcell; 2 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.

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.
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 domelike 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 Dz. 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 Dz. 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 a, while the last phase (3 -1) is related to a decreasing a.

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, 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: threedimensional 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).
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 o t o 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 o t o 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 4 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 Dz CoM = À6.6 Â 10 À4 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 p À(Dz/R out ) 3 for small prescribed shape heights Dz and collapses onto the results for the swelling disk.
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. 50 We can expand the axisymmetric flow field in the far-field into Legendre polynomials 11,12 and extract the dipole contribution p 2 P 2 (cos y)/r 2 of the radial part u r (r,y) of the flow field u( r) from eqn (23) and normalize by the center of mass velocity |v CoM | to define the dimensionless parameter g 2p 2 3R out 2 jv CoM j : Values g o 0 (p 2 o 0) indicate ''pulling'' motion of the swimmer, values g 4 0 (p 2 o 0) ''pushing'' motion, and g E 0 (p 2 E 0) a neutral motion. 3 The results for g 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 shortwavelength 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.

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 a and the outer annulus with 1/a, we saw a transition into an elliptic dome-like shape for a 4 1 and a transition into a hyperbolic saddle shape for a o 1. The control parameter of this shape transition is the swelling factor a. 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 a 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 a. 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. 2 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 (|Dz CoM /R out | B À(Dz/R out ) 3 ). 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. 31 These disks had a diameter of 30 mm and a thickness of 5 mm. 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 mm and a swimming distance of about 1% of the radius after ten cycles, a net velocity of 0.075 mm s À1 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 (|Dz CoM | p h 3 ).
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 t h B ZR out /s 0 B g FvK 3/2 Z/Y 3D , which is of the order of t h B 10 À5 s for the thermoresponsive hydrogel disks from ref. 31. Then Re o 1 is equivalent to t h 4 rR out 2 /Z and restricts the low Reynolds number regime to sufficiently thin disks with h 3 /R out o 10 À1 mm 2 . 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 51 or shell-based snapping systems. 52 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. 26 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.