Method of actuating a shape changeable member, shape changeable member and actuating system The present invention relates to a method of actuating a shape changeable member of actuatable material. The invention further relates to a shape changeable member and to a system comprising such a shape changeable member and a magnetic field apparatus. Miniature members are envisaged that can facilitate access to and inspection of narrow, remote and unexplored environments. However, restricted mobility of miniature members may hinder the realization of targeted applications in highly unstructured landscapes. To date the miniature members produced are generally tethered members, this means that they are generally connected to a further structure and can be adapted to move relative to the further structure. While there also exists untethered miniature members, these miniature members are generally only capable of carrying out one mode of locomotion and are specifically designed to carry out that one single mode of locomo- tion. The use of members that can only carry out limited number of modes of locomotion however hinders the movement capabilities and applications that can be carried out by such members.

In view of the foregoing it is an object of the present invention to make available shape changeable members that can more readily access remote and unexplored environments and that are capable of carrying out multiple modes of locomotion.

This object is satisfied in accordance with a method of actuating a shape changeable member as defined in claim 1 .

By way of example, such a method of actuating a shape changeable member, in which the shape changeable member has a height, a width as well as a length, and is untethered; comprises the steps of:

- actuating the shape changeable member by means of at least one stimulus to bring about a change of shape of the member at least over the length of the shape changeable member in at least one spatial direction and thereby inducing a torque into the shape changeable member; - effecting a locomotion of the member by means of the induced change of shape of the member in the at least one spatial direction; and

- moving the shape changeable member in the at least one spatial direction. The use of untethered members increases the movement capabilities of such members. Moreover, the use of at least one stimulus to bring about a variety of changes of shape of the member to induce different moments into the member to effect different types of locomotion increases the variety of modes of moving the member. The different modes of moving the member translate to different types of locomotion which were previously not possible using one and the same member.

Thereby miniature mobile members capable of accessing and navigating effectively in narrow remote environments can be made available. These miniature members may enable unprecedented applications in microfactories, bioengineering, environmental remediation, and healthcare. Moreover, one and the same miniature device can be made available that can inherently carry out the different modes of locomotion.

In this connection it should be noted that the moment can also be referred to as torque and the induced torque is generally small when the induced change of shape of the member causes a small deformation of the member. Moreover, it should be noted that the torque vanishes when the net magnetic moment of the deformed member is aligned with an applied stimulus, such as an applied magnetic field.

Preferably the induced torque is used to steer the shape changeable member in the at least one direction.

It is preferred if a type of locomotion of the shape changeable member is selected from the group of members consisting of swimming, walking, flipping, crawling, rolling, diving, immersion, emersion, jumping, landing, surface climbing, liquid meniscus climb- ing, hopping, crawling inside a confined space, spinning, flying, gliding.

Dexterous miniature members that can swim, walk, roll, jump, or crawl have been proposed in recent years. However, these members are only capable of performing one specific locomotion mode, and are inadequate to escape from liquid environments to solid terrains. In contrast to this the shape changeable member can be actuated to carry out several modes of locomotion one after the other. This means that the member can, for example, move through a liquid, leave the liquid and then move on a solid surface, and vice versa. For example, untethered, amphibious, millimeter-scale soft (rubber like) members capable of multiple modes of locomotion can be made available that move across and/or through liquid and solid unstructured terrains. In this connection, soft means that the member is inherently flexible, such as a member comprising rubber and is not a hard/rigid member such as a piece of non-flexible metal or wood.

The members can be actuated by locomotion mode-specific external magnetic fields, the magneto-elastic sheet-shaped members with a pre-programmed magnetization profile and hydrophobic surface can swim inside and on the surface of liquids, climb liquid menisci, roll, walk on ground, jump directionally over obstacles, and crawl within narrow spaces. Such members can reversibly transit from the surface to the bulk of a liquid, and from a liquid surface to rigid ground.

Advantageously the change of shape comprises at least one of a deformation, a contraction, a flexing, an undulation and a stretching of the actuatable material of the shape changeable member. In this way, the member can change its shape in a variety of ways with each change of shape being capable of being used for certain types of locomotion.

By way of example, the shape changeable member can be deformed to adopt a sinelike or cosine-like shape to affect a walking type of locomotion, in particular wherein a rotating magnetic field B sequence of small magnitude is applied to the shape changeable member to enable an undulating locomotion of the shape changeable member, with an applied magnetic field preferably having a small magnitude |B| being selected in the range of 0.1 to 200 mT in particular of 2 to 5 mT and especially of 3 mT. By way of a further example, the shape changeable member can be deformed to adopt a part-circular shape, such as a semi-circular shape, to effect a rolling or meniscus climbing type of locomotion, in particular wherein a rotating magnetic field B sequence of large magnitude is applied to the shape changeable member to enable a rolling locomotion of the shape changeable member, with the applied magnetic field preferably effecting a rolling type of locomotion in a clockwise sense, and with an applied magnetic field preferably having a large magnitude |B| selected in the range of 0.1 to 1000 mT, in particular of 12 to 20 mT and especially of 15 mT, in particular wherein, to effect a meniscus climbing type of locomotion a large constant magnetic field B is applied followed by a rotating magnetic field B sequence of large magnitude with an applied magnetic field preferably having a large magnitude |B| being selected in the range of 0.1 to 1000 mT, in particular of 10 to 20 mT and especially of 12 mT.

By way of a further example, the shape changeable member can be deformed to un- dulate and generate a travelling wave on the member to affect a swimming type of locomotion or crawling type of locomotion.

By way of a further example, the shape of the shape changeable member can be flipped from a first shape into a second shape and back into the first shape to affect a jelly-fish swimming type of locomotion, wherein the first and second shapes preferably resemble that of a part-circle.

By way of yet a further example, the shape changeable member can be pulsed with the at least one stimulus in order to affect the jumping type of locomotion.

By way of yet a further example to affect an immersion type of locomotion a magnetic field having a rotating B sequence is applied, with the magnetic field effecting a rotation of the shape changeable member in a counter-clockwise sense to reduce a contact with a liquid surface, with the applied magnetic field preferably having a magnitude selected in the range of 1 to 200 mT, in particular of 15 to 30 mT, especially of 20mT.

By way of yet a further example to affect a landing type of locomotion a magnetic field having a rotating B sequence of variable magnitude applied, with the initial B- sequence having a high magnitude which is then reduced to a low magnitude, with the high magnetic field preferably being selected in the range of 1 to 1000 mT, especially of 10 to 15 mT and the low magnitude being selected in the range of 0.1 to 100 mT, especially of 3 to 8 mT.

By way of yet a further example a rocking sequence of the magnetic field is applied to perform the walking type of locomotion, in particular the walking type of locomotion is effected by use of four phases in which four different kinds of magnetic field are applied.

By way of yet a further example to affect a jelly fish-like swimming motion mode of lo- comotion of the shape changeable member the applied magnetic field is sequentially flipped from a positive magnetic field to a negative magnetic field causing an inversion of the shape of the shape changeable member.

Preferably the stimulus is one of an external stimulus and an internal stimulus. Examples of an internal stimulus can e.g. be pressurized air that is provided at the member in order to inflate or deflate certain regions of the member to bring about an actuation of its shape. Examples of an external stimulus are the application of a magnetic field, an electric field, heat, or laser light. E.g. direct illumination of the laser on the member which causes ejection of material and motion of the member.

Examples of members that can be actuated by external stimuli that do not include the application of a magnetic field as the main source of stimulus are thermally swelling materials, light-sensible materials locally heated by laser light. For example, shape- memory polymers can be used that can deform when they are stimulated by heat. By pre-programming these polymers, different regions of the polymer can create different curvatures.

It is however preferred if the external stimulus is present in the form of a time varying magnetic field applied in at least one dimension, preferably two dimensions, most preferably in three dimensions. A time varying magnetic field is simple to realize and can be implemented in a cost effective manner and can be used to move a variety of members through various types of media.

Moreover, a time varying magnetic field can not only be changed in its magnitude but also in its direction and this variation in field strength and direction can be utilized by a shape changeable magnetic member that has a pre-determined magnetization profile to change its shape over its length and height in quick succession in order to be moved by the magnetic field using different forms of locomotion. Preferably a magnitude of the time varying magnetic field applied is selected in the range of 0 to 1000 mT, especially of 0 to 500 mT. Such magnetic fields can be used to move a member both in liquid environments as well as in dry environments.

Advantageously the magnitude of the time varying magnetic field applied is varied with a frequency selected in a suitable range of, e.g., 0 to 1000 kHz, preferably of 0 to 200 Hz, especially of 0 to 50 Hz. These frequencies have found to be advantageous to bring about the locomotion of the shape changeable member.

It is preferred if a direction of the time varying magnetic field applied is varied with a frequency selected in the range of 0 to 100 kHz, especially of 0 to 100 kHz. These fre- quencies have been found to be advantageous to bring about the change of shape of the member and to induce the various types of locomotion.

It is preferred if the shape changeable member can be actuated in order to additionally carry out at least one function. In this way, for example, a robot can be made available that is capable of executing pick-and-place tasks and shape deformation-triggered release of cargos.

For example, the at least one function comprises at least one of a gripping function, a holding function, a dropping function, a clamping function and picking up function.

In accordance with a further aspect the invention also relates to a shape changeable member of actuatable material, the shape changeable member being untethered and having a width, a height and a length, wherein the shape changeable member is con- figured to be actuated by means of at least one stimulus in order to bring about a change of shape of the member at least over the length of the shape changeable member; wherein the shape changeable member is configured to be actuated in accordance with a method of actuating the shape changeable member as discussed in the foregoing. In this connection it should be noted that the shape changeable member can be furthered using similar features to those discussed in the foregoing in relation to the method in accordance with the invention.

Such minimalist, programmable and versatile soft members could enable applications in environmental monitoring and minimally invasive medical operations inside the hu- man body.

The proposed members with versatile and robust multimodal locomotion capability may be tailored to diverse applications, including environmental monitoring and minimally invasive medical operations. Further progress is envisaged by implementing a remote member localization method and improving the member's dexterity, precision, safety, and functionality using a closed-loop feedback control by monitoring the member while this is moved e.g. through the human body.

Preferably the stimulus is an external stimulus in the form of a magnetic field applied in at least one dimension and most preferably in three dimensions, with the magnetic field interacting with a pre-programmed magnetization profile of the shape changeable member, with the member preferably having a hydrophobic surface. Providing a magnetized member that moves in a magnetic field is a beneficial way of implementing the various modes of locomotion required to e.g. facilitate drug delivery in the human body. A viable solution to develop miniature devices with multimodal and multi-terrain locomotion capability is to employ soft active materials, such as ecoflex, shape memory alloys, liquid metals, silicone rubbers, silicone-based materials, polyurethanes, soft gels (hydrogels, oil-based gels, aerogels) as well as natural polymers, such as DNA or proteins that are embedded with ferromagnetic, paramagnetic or diamagnetic parti- cles, which endow machines with extra degrees-of-freedom and support programmable shape changes while assuring safe interaction with delicate environments.

Moreover, magnetic actuating fields can harmlessly penetrate most biological materials, making them particularly promising for biomedical applications.

In a further aspect the present invention relates to an actuating system comprising a shape changeable member in accordance with the foregoing and a magnetic field apparatus, wherein the magnetic field apparatus is a source of the at least one stimulus. The advantages discussed in the foregoing in connection with the method of actuating as well as in connection with the shape changeable member likewise apply to the system.

It is preferred if the magnetic field apparatus is configured to generate a uniform 3D time-varying magnetic field at a center of the system. In this way the shape changeable member can be actuated on by a magnetic field in order to effect the different types of locomotion. Such as system can then be used to e.g. steer miniature robots to different positions of application. Preferably the magnetic field apparatus comprises at least two, preferably three and especially four, pairs of magnetic coils, wherein each of the pairs of coils is arranged on a common axis, with the respective common axes being arranged at least approximately perpendicular to one another at a center of the system. Such a coil arrangement can generate very accurate time-varying magnetic fields that are suitable to moving the shape changeable member.

It should be noted in this connection that if four pairs of coils are used then the fourth pair of coils is configured to control a gradient of the applied magnetic field. Advantageously the shape changeable member is configured to be actuated by means of the pairs of magnetic coils in such a way that the shape changeable member carries out at least one type of locomotion within the applied magnetic field. It is preferred if each coil of each pair of magnetic coils is arranged at least substantially equidistant from the center. In this way a homogenous magnetic field can be generated that can be varied in a time dependent manner both with regard to the direction of the magnetic field applied and with the magnitude of the applied magnetic field. In this way the pairs of magnetic coils can be configured to produce a time-varying magnetic field at the center.

In order to produce a time-varying magnetic field at the center, each pair of coils can be actuated to produce a magnetic field at the center whose field strength can be suit- ably varied, preferably in the range of 0 to 500 mT, in particular of 0 to 50 mT. Indeed each pair of coils can be actuated to produce a magnetic field at the center whose field strength can be varied with a frequency in the range of at least 0 to 1 kHz, in particular of 0 to 50 kHz. In fact each pair of coils can be actuated to produce a magnetic field having a varying direction in one of the following planes spatial planes X-Y, X-Z, Y-Z. By changing the direction in the different planes in which the field is applied means that the shape changeable member can be moved in all three directions of space. Preferably the change in direction of the magnetic field applied can be varied with a frequency selected in the range of 0 to 100 kHz, in particular of 0 to 200 Hz.

Advantageously each coil of each pair of coils is configured to carry a current in the range from 0 to 20 A for a coil comprising 3 layers of wires with 60 wires per layer and a diameter of 1 .05 cm.

Preferably, if three pairs of magnetic coils are used, two pairs of the magnetic coils are provided with an at least substantially equal size and each coil of the third pair of magnetic coils has a larger size than a respective coil of the other two pairs of magnetic coils.

In this way, a simple configuration is made available to produce a uniform magnetic field. Advantageously a diameter of each of the four coils of the two pairs of magnetic coils of at least substantially equal size is selected in the range of 50 to 65 cm and a diameter of each of the two coils of the third pair of magnetic coils of at least substantially equal size has a size that is selected in the range of 100 to 140 cm.

Preferably the shape changeable member is configured to be steered by a magnetic field applied in at least one of the Y-direction and the Z-direction and is further configured to be changed in shape by a magnetic field applied in the X- and Y- direction and is configured to be moved by the magnetic fields applied in the X- and Y- direction.

Advantageously the magnetic field is applied relative to a local X-Y frame of the shape changeable member and through a presence of the magnetic field deforms a shape of the shape changeable member in the X-Y frame, wherein the shape changeable member is steering by a magnetic field applied in at least one of the Y-direction and the Z-direction.

The invention will be explained in detail in the following by means of embodiments and with reference to the drawing in which is shown: Fig. 1 a) a design of a magneto elastic-shape changeable actuatable member and b) the magnetic response of the magneto elastic-shape changeable actuatable member;

Fig. 2 various types of modes of locomotion and associated transition modes of the magneto elastic-shape changeable actuatable member, with a) illus- trating swimming, b) meniscus climbing, c) immersion, d) landing, e) rolling, f) walking, g) crawling and h) jumping;

Fig. 3 multimodal locomotion over liquid-solid hybrid terrain by the magneto elastic-shape changeable actuatable member, with a) illustrating rolling, diving and swimming, b) immersion, c) emersion, d) meniscus climbing, landing, jumping, walking and e) walking and crawling;

Fig. 4 a) approach and pick up of a load by the magneto elastic-shape changeable actuatable member and b) transport and placement of the load by the magneto elastic-shape changeable actuatable member;

Fig. 5 a magnetization process of a soft magneto-elastic shape changeable actuatable member, with a) illustrating an elastomeric beam loaded with

NdFeB microparticles wrapped around a cylindrical glass tube (perimeter: 3.7 mm) and magnetized by a uniform 1 .5 T magnetic field produced by a vibrating sample magnetometer, and b) illustrating the programmed magnetization profile in the unwrapped member, for _R = 45°; Fig. 6 a six-coil electromagnetic setup to actuate the magneto-elastic shape changeable actuatable member;

Fig. 7 a quasi-static analysis of a soft shape changeable actuatable member;

Fig. 8 a rotating B sequence enabling the undulating motion of the magneto- elastic shape changeable actuatable member, with (a) illustrating the required B in dependence on time to undulate the member 10 to the right, with the variables B _x and B _y representing the magnetic field along the x- and y-axis of the coordinate system, shown in (b) with Fig. 8(b) illustrating two principle shapes;

Fig. 9 further illustrations as shown in Fig. 8a) and b) for a further type of locomotion;

Fig. 10 further illustrations as shown in Fig. 8a) and b) for a further type of locomotion;

Fig. 1 1 further illustrations as shown in Fig. 8a) and b) for a further type of loco- motion;

Fig. 12 further illustrations as shown in Fig. 8a) and b) for a further type of locomotion;

Fig. 13 further illustrations as shown in Fig. 8a) and b) for a further type of locomotion;

Fig. 14 further illustrations as shown in Fig. 8a) and b) for a further type of locomotion;

Fig. 15 snapshots of straight and directional jumping by the magneto-elastic

shape changeable actuatable member, with a) illustrating a first type of straight jumping, b) a second type of straight jumping, c) a type of direc- tional jumping and d) the magnetization profile of the shape changeable magnetic member used in a to c, and its magnetic response at large |B|; Fig. 16 a vibrational analysis of a shape changeable actuatable member; Fig. 17 boundary conditions for the Type-1 straight jump of Fig. 15a;

Fig. 18 a design of a modified shape changeable actuatable member for a mag- netically triggered cargo release, with a) illustrating an extra strap added to the member 10, b) the extra strap locked into a hole on the body of the member with cargo arranged there between and c) the opening of the strap;

Fig. 19. (a) The rectangular sheet-shaped magnetic soft robot's dimensions and magnetization profile, m. (b) Theoretically predicted and experimental shapes of a magnetically deformed millirobot. a is a clockwise angle from the x-axis, and it is used for describing the direction of B _xy . (I) The robot is in its rest state under null B _xy . A small residual curvature is induced by the pre-stress resultant from the demolding process. (Il-lll) Small defor- mation shapes under a small magnitude B _xy {B _xy = 5 mT) aligned along the two principal directions (a = 225° for sine and 315° for cosine shape) in the xy plane. In contrast, the robot deforms into a 'C'-shape (IV) or 'V- shape (V) when a large magnitude B _xy (B _xy = 20 mT) aligns along the principal axis (o = 135° and 315°, respectively). The induced pre-stress from the demolding process may introduce a deviation between the predicted and actual shapes, particularly for (II) and (V). The theoretical M _ne t of the robot for each shape is shown on the right, (c) The clockwise rotation of the 'C'-shaped robot is induced by rotating B _xy . Scale bars: 1 mm. Fig. 20 the sub-panels for the swimming locomotion in (a) are displayed in a horizontal sequence to show the relative vertical displacement between them. The sub-panels for the rest of the locomotion are arranged in a vertical sequence to show the relative horizontal displacement between them. The corresponding time stamps and B in each sub-panel are shown in the bottom-right and top-left corners, respectively, (a) Jellyfishlike swimming in water using a time-symmetrical but speed-asymmetrical gait, (b) Water meniscus climbing. A counter-clockwise magnetic torque progressively adapts the pose of the robot as it deforms and ascends due to buoyancy, (c) Landing, i.e., transition from water surface onto solid ground. A clockwise rotating B peels the robot off the water surface and lets it stand on the platform, (d) Immersion, i.e., transition from the surface into the bulk of a water pool by a combination of curling and rigid- body rotation, (e) Rolling by a clockwise rotating B of high magnitude, (f) Walking. The robot tilts and changes its curvature to create a net stride in each cycle, (g) Crawling inside a tubular channel with a cross section of 0.645 x 2.55 mm ² by using an undulating traveling wave along the robot's body, (h) Directional jumping. The robot uses its rigid-body motion and shape change to induce a jumping momentum. More than one robot is used for this illustration but we highlight that all these robots have the same design (scale bars: 1 mm); and

Fig. 21 (a) the soft robot navigating across a synthetic stomach phantom using a combination of locomotion modes, (b) Ultrasound-guided locomotion: the robot (marked out by dotted lines) rolls within the concealed areas of an ex vivo chicken tissue, as visualized by ultrasound imaging, (c) The robot approaches a cargo item (nylon, 1 χ0.8*1 .5 mm ³ ) by walking on a flat rigid surface, picks up the cargo by curling into the 'C'-shape, transports away the cargo by rolling and maintaining its 'C'-shape configuration, and releases the cargo by uncurling at a new position, (d) Dynamic and selective cargo release. A paper tissue (0.5x0.5x0.1 mm ³ , used as model drug container) is bound to the robot by an extra appendage After pre- bending the robot, B is quickly reversed to open the appendage and re- lease the cargo (cargo highlighted by the white dashed line) (scale bars: 1 mm).

Fig. 1 a shows a design of a shape changeable actuatable member 10, present in the form of a strip of material . The shape changeable actuatable member 10 can also be referred to as a magneto-elastic member 10. The member 10 has a length I in the x- direction, a height h in the y-direction and a width w in the z-direction. In the present example, the member 10 is formed from a continuum of infinitesimal segments ds of which nine segments s are shown for the purpose of illustration along the x-direction, with the segments s (ds) each comprising the same material but having a different direction of magnetization as indicated by arrows.

The final segment Sf in the x-direction indicates a phase shift β _κ of the magnetization profile, this phase shift is important as will be discussed in the following.

The member 10 is magneto-elastic and fabricated from a hydrophobic silicone rubber sheet loaded with embedded 5 μηη NdFeB particles. Fig. 5 and the subsequent description details steps that can be carried out to form such a member 10. The member 10 is magnetized as a single-period sinusoidal function for its m _x and m _y component along the sheet's length I and equally along the sheet's width w and thickness h.

Fig. 1 b shows the member 10 of Fig. 1 in four different states of excitation. The different states of excitation are brought about by way of remote magnetic actuation by means of a multi-coil electromagnet system (see Fig. 6 in this regard). This multi-coil electromagnet system enables untethered operations, i.e. the member 10 is not connected to a host, substrate or the like and can indeed move freely in space.

The first state shows a member 10 with β _κ = 45 ° in its rest state when no magnetic field B is applied, a scale bar is also indicated and has the size of 1 mm. The member 10 shown in the present example has a size of 3.7 χ 1 .5 χ 0.24 mm ³ .

Typical sizes of the member can be selected in the range of 0.001 to 100 χ 0.001 to 100 x 0.001 to 20 mm ³ especially of 0.01 to 10 χ 0.01 to 10 χ 0.01 to 2 mm ³ . It should be noted in this regard that also non-uniform shaped members 10 can be produced; however, for the ease of simplicity a member 10 of strip shape is described in the following. The second state shows an applied magnetic field of 3 mT field strength that is aligned in the X-Y plane with a principal direction of a = 225°. In the second state the magnetic member 10 adopts a sinusoidal shape. The third state shows an applied magnetic field of 3 mT field strength that is aligned in the X-Y plane with a principal direction of a = 135°. In the third state the magnetic member 10 adopts a cosine like shape, i.e. by changing the phase of alignment of the magnetic field by 90° has caused a shift in the general shape of the member by 90°. If one were to continue changing the principal direction of the magnetic field of the same amplitude one can observe a traveling wave traversing through the member 10. In this way, the member 10 can be caused to undulate and thereby a swimming like type of locomotion can be effected (see also Fig. 3a) In this regard, it should be noted that a time varying magnetic field can be applied with a field strength selected in the range of 0 to 1000 mT, especially of 0 to 500 mT, with a frequency of the magnitude being variable in a range selected as 0 to 1000 kHz, especially of 0 to 1 kHz. Likewise, a direction of the applied magnetic field can be changed with a frequency selected in the range of 0 to 1 kHz.

The fourth state depicted shows the application of a significantly larger magnetic field in comparison to the third state for the same principal direction of magnetization. In this magnetic field, the member 10 is distorted into a horse shoe like shape or a semicircle like shape when |B| is large.

The four states show that the shape of the member 10 can be changed by changing the principle direction and/or the strength of the applied magnetic field. Indeed, by generating desired magnetic fields net magnetic moments arise in the member 10. The electromagnets allow the soft member 10 to thus not only produce time-varying shapes, but also various magnetic moments that can invoke a variety of modes of locomotion. The magnetic moments align with an externally applied magnetic field B, so that magnetic torques can steer the member 10 in three-dimensions (3D) in all modes of locomotion. Fig. 2 shows various types of modes of locomotion and associated transition modes of the magneto elastic-shape changeable actuatable member 10. The scale bars shown in Fig. 2 relate to an actual size of 1 mm.

Fig. 2a shows a mode of locomotion of the member 10 that resembles Jellyfish-like swimming, when the member 10 is completely immersed in water. Indeed, it has been found that the member 10 can overcome gravity by swimming upwards. The swimming is brought about by a time-asymmetrical reversal of B along the vertical direction. The maximum field strength used is |B| _max = 30 mT. Thus, a periodic magnetic field B with a time-varying magnitude along the vertical direction causes the jellyfish-like swimming motion of the member 10, and propels the member 10 to reach the water surface (Fig. 3c). Upon emersion, the soft member 10 strongly pins to the water-air interface by exposing its hydrophobic body surface to air.

Fig. 2b shows a mode of locomotion of the member 10 that resembles water meniscus climbing. A counter-clockwise magnetic torque progressively adapts the pose of the member 10 to the meniscus profile and raises its center of mass. The maximum magnetic field used in this example is |B| _max : 12 mT.

By applying the different magnetic field a programmed shape change occurs from the shape shown in Fig. 2a that allows the member 10 to exploit the capillary action of a water surface and move upward along the slope of a concave meniscus (Fig. 2b). Meniscus climbing is typically adopted by insects when they cannot swim across friction- less liquid barriers, thus, the member 10 can adopt an insect like mode of locomotion. It has namely been found that a control of the magnitude of B and of the direction of B makes the member 10 curve upwards and gradually rotate to match the local curvature of the water surface. This curving and rotation permits the member 10 to undergo capillary displacement without extra energy expenditure. Conversely, a fast sequence of downward bending and rotation disengages the member 10 from the water-air inter- face as shown in Fig. 2c.

Fig. 2c shows a mode of locomotion of the member 10 that resembles immersion, i.e., a transition from the surface of a water pool into the bulk water of a water pool by a combination of curling and rigid-body rotation. The maximum magnetic field used in this example is |B| _max : 30 mT.

Fig. 2d shows a mode of locomotion of the member 10 that resembles landing, i.e., a transition from a water surface onto solid ground. A magnetic field that changes clockwise B peels the member 10 off the water surface and lets it stand on the platform. The maximum magnetic field used in this example is |B| _max : 12 mT. On meniscus climbing, the member 10 can thus efficiently rotate to peel away from the water surface to stand on an adjacent solid substrate. The low wettability (receding water contact angle -78°) of the microscopically rough elastomer surface of the member 10 leaves the receding triple contact lines unpinned and makes water easily dewet the surface of the member 10 during its displacement. Fig. 2e shows a mode of locomotion of the member 10 that resembles rolling by means of a high magnitude, and low frequency rotating B field. The maximum magnetic field used in this example is |B| _max : 20 mT. By rolling, i.e. rigid-body rotation around the z-axis enabled by means of the high magnitude, low frequency B, the member 10 can advance directionally over a rigid substrate and also dive from a solid onto a liquid surface (see Fig. 3a). However, the curled member 10 cannot roll across substrate gaps wider than the diameter of the member 10. In order to bridge gaps the member 10 changes into a walking mode of locomotion.

Fig. 2f shows a mode of locomotion of the member 10 that resembles walking. The direction and magnitude of the magnetic field B control the gait of the member 10 and modulate the stride between its legs, respectively. The maximum magnetic field used in this example is |B| _max : 7.7 mT. The dashed lines mark the net displacement after one walking cycle. Walking is a particularly robust way to move over unstructured terrain morphologies, and affords precise tuning of stride granularity. Directional walking can be achieved by sequentially and periodically adapting the tilting angle and curvature of the member 10 through the control of the orientation and magnitude of the magnetic field B, respectively.

By way of example the walking to the right shown in Fig. 2f will be described in more detail by using the following sequence of the magnetic field B (see also Fig. 13a):

1 . Start walking (Frame I): The magnetic field B is oriented at a=135° with a magnetic field strength |B| = 7.7 mT. Under this condition, the member 10 assumes a vertical- ly symmetric pose with both feet touching the ground (see also Fig. 12b).

2. Rise of the front foot (Frames I and II): B rotated from a = 135° to 124° while |B| = 7.7 mT. This lifts up the front foot (with respect to the moving direction) while keeping the back foot fixed.

3. Stretch and lowering of the front foot (Frames II and III): While a is kept at 124°, |B| is linearly decreased from 7.7 mT to 0 mT. This lowers the front foot and extends the horizontal span of the member 10 by an amount d while the member 10 pivots on the back foot.

4. Shrink and rising of the back foot (Frames III and IV): While a is kept at 146°, |B| increases from 0 mT to 7.7 mT. This makes the member 10 curl to shrink its span and lift up the back foot while keeping the front foot fixed as a pivot.

5. Lowering of the back foot (Frames IV and V) (Lower the back foot): B is rotated from a =146° to 135° while |B| = 7.7 mT. This closes the control cycle by returning the member 10 to its starting pose at a net rightward displacement d with respect to its initial position (compare to frame I). Jellyfish-like swimming (Fig. 13)

The jellyfish-like swimming shown in Fig. 2a is realized by the following steps:

1 . The member is acted on to implement a fast power stroke.

2. The member is then acted on to implement a slow recovery stroke.

By this periodic flipping of B along the vertical direction to give a fast power (downward) stroke, the robot (member), functioning at a Reynolds number Re -32.5, is shown to swim upward against the gravitational field, and even against the downward magnetic field gradient eventually present at the vicinity of the bottom coil of the work- ing space.

Meniscus climbing (Fig. 9)

The meniscus climbing shown in Fig. 2b is realized by the following steps:

1 . The member is acted on to assume a vertically symmetric pose with upward curva- ture. The member achieves in 400 ms a stable position along the meniscus known to depend on both the initial meniscus profile and the curvature of the device.

2. B is subsequently rotated clockwise toward to progressively displace the position of the robot along the meniscus and raise its center of mass. Once the device reaches contact with the solid substrate, B is turned off and the robot rests on the edge of the adjacent platform (Frame I in Fig. 2d).

Immersion (Fig. 10)

The immersion shown in Fig. 2c is realized by the following steps:

1 . The member is acted on to initially curl downward (Frame II in Fig. 2c).

2. The member is then acted on to rotate counterclockwise to reduce the contact with the water-air interface to a minimum (Frame II in Fig. 2c). At the end of this step, one ending of the beam is still attached to the water surface (Frame III in Fig. 2c).

3. A quick 180° flipping of the external B field disengages the ending of the member from the water surface and let the robot sink (Frame IV and V in Fig 2c).

Landing (Fig. 1 1 )

The landing shown in Fig. 2d is realized by the following steps:

1 . The member is acted on to keep rotating clockwise while using the edge of the platform as the pivot point.

2. When the second leg of the member reaches the platform (Frame V, Fig. 2d), the external B field is turned off.

Rolling (Fig. 8) The rolling shown in Fig. 2e is realized by using the following steps:

1 . A large magnitude B (15 mT) along a field angle of 157° to deflect the robot into a distorted cosine shape resembling a semicircle-like shape (Fig. S5b).

2. The extermal B starts rotating, a rigid-body magnetic torque, Xrigid-body, is induced by the robot, which can be described mathematically as:

The external B starts rotating and the agent, the robot's M _ne t will be aligned with its final direction, and will be reduced to a null vector.

Walking (Fig. 12)

The rightwards walking shown in Fig. 2f is realized by the following steps

1. Start of walk (Frame I): B oriented at a=135° with |B| = 7.7 mT. Under this condition, the robot assumes a vertically symmetric pose with both feet touching the ground (Fig. S8b).

2. Rise the front foot (Frame I and II): B rotated from a = 135° to 124° while |B| = 7.7 mT. This lifts up the front foot (with respect to the moving direction) while keeping the back foot fixed.

3. Stretch and Lower front foot (Frame II and III): While a is kept at 124°, |B| is linearly decreased from 7.7 mT to 0 mT. This lowers the front foot and extends the horizontal span of the robot by an amount d while the robot pivots on the back foot.

4. Shrink and Rise the back foot (Frame III and IV): While a is kept at 146°, |B| increases from 0 mT to 7.7 mT. This makes the robot curl to shrink its span and lift up the back foot while keeping the front foot fixed as pivot.

5. Lower the back foot (Frame IV and V) (Lower the back foot): B is rotated from a =146° to 135° while |B| = 7.7 mT. This closes the control cycle by returning the robot to its starting pose at a net rightward displacement d with respect to its initial position (compare to frame I). Crawling and undulating swimming (Fig. 7)

The crawling mode of locomotion (Fig. 2g) require undulating, time-varying shapes. The agent achieves its undulation by using a rotating B of low magnitude and relatively high frequency. The deformation produced by the time-dependent undulation can be described mathematically as: where the variable T represents the period of the rotating B. The same motion is also applied to undulating swimming.

Jumping

The jumping shown in Fig. 2h is realized by using the following steps:

1 . In Type-1 straight jumping, the agent is first pre-bent by the magnetic field to curl upward and stand on its center (frame I of Fig. 2h).

2. The direction of B is then quickly reversed to force the agent to flatten out (frames II and III of Fig. 2h) and hit the substrate by its feet.

Fig. 2g shows a mode of locomotion of the member 10 that resembles crawling inside a cylindrical tube (inner diameter: 1 .62 mm). A clockwise rotating magnetic field B (B = 3 mT at 15 Hz) propels the motion to the right by body undulation. The dashed lines mark net displacement after five crawling cycles. When facing narrow openings in standing obstacles that may prevent walking, the member 10 can effectively resort to small amplitude body undulations to crawl through the opening. The crawling motion is encoded by a low magnitude rotating magnetic field B that produces reversible longitudinal traveling waves along the magneto-elastic member 10. A similar control se- quence additionally makes the member 10 swim efficiently on liquid surfaces akin to a Taylor swimming sheet (see Fig. 3a), as it does in the member 10 with m featuring multiple sinusoidal periods.

Fig. 2h shows a mode of locomotion of the member 10 that resembles directional jumping. The member 10 is initially tilted at 9° by the applied magnetic field B. The magnetic field B is then flipped clockwise and maintained for 25 ms. The maximum magnetic field used in this example is |B| _max : 18.9 mT. The arrows indicate the orientation of B in the XY-plane and the instantaneous net magnetization vector of the member 10, respectively.

Fronted by obstacles too high to roll or walk over, the soft member 10 can thus also resort to jumping by imparting an impulsive impact on a rigid substrate. Both straight jumping (see Fig. 15) and directional jumping strategies (see Fig. 2h and Fig. 15) can be achieved. To steer the jumping mode of locomotion, the control sequence of the magnetic field B prompts both the rigid-body rotation of the member 10, which specifies the jump direction, and elastic deformations to maximize the momentum of the member 10 before striking the substrate. Alongside the geometrical member design and the control parameters of the magnetic field B, the net magnetic moment of the deformed member additionally allows the modulation of jump height by magnetic gra- dient pulling. Magnetic gradient pulling alone could theoretically levitate the member 10; however, the required gradient magnitude is too large to be made practically feasible.

Fig. 3 presents a hybrid liquid-solid environment that a member 10 can fully explore only through the combination of all of its modes of locomotion. A long-ranged water meniscus adjacent to a smooth, hydrophilic solid platform prompts capillary climbing by the member 10, which subsequently jumps directionally beyond an otherwise insurmountable, fixed obstacle. Navigation across the unstructured environment of e.g. a surgical human stomach phantom further illustrates the robustness of the multimodal member locomotion towards future potential biomedical applications.

The environment shown in Fig. 3 features laser-milled, poly(methyl methacrylate) platforms (water contact angle < 10°) in a Plexiglas box (52 χ 32 χ 25 mm ³ ) containing de- ionized water. The various modes of locomotion shown in Fig. 3 are sequentially cap- tured in separate videos of which frames are shown, the scale bars used in these Figs, scale to 1 mm.

In the following the procedures to implement the multimodal mode of locomotion of the member 10 shown in Fig. 3 will be described. In Fig. 3a, the member 10 first uses its rolling locomotion to get onto the water surface. After dipping into the water, the member 10 drifts slightly away from the platform along the meniscus and then starts swimming after it has stabilized. A clockwise magnetic field B (|B| = 3 mT, 25 Hz) is applied to let the member 10 swim to the right. Subsequently, the member 10 sequentially bends downwards on the water surface to disengage from the water surface. The member 10 is then and rotated counterclockwise to sink into the water, and subsequently swims back to the water surface to engage again at the water/air interface (Fig. 3b-c).

Fig. 3d shows the member 10 using four different modes of locomotion on land. After climbing onto the solid substrate by way of meniscus climbing the member 10 is present on the solid platform. The net magnetization of the member 10 is subsequently used to rotate the member 10 to the required initial pose for directional jumping to move the member 10 beyond the otherwise unsurmountable obstacle. After this, the member 10 starts to walk.

It has been observed that the member 10 incurs slight slipping before and after jumping across the obstacle due to the magnetic gradient force. As in this case, slip may happen because the currently implemented control loop is open: The member 10 traverses a large working space, and for each position in the space the suitable mag- netic actuation matrix is chosen by the human controller. Thus, a small gradient force associated with a control error may be present. Specifically, the slipping happens when the member 10 has only one foot in contact with the substrate, that is, when the friction is smallest and even a small gradient force can displace the member 10. It is worth noting that the walking can make the member 10 advance against an opposing magnetic gradient pulling force. Moreover, as the net magnetization of the member 10 can steer the member 10, it can be seen that the walking member 10 is initially misaligned with the desired advancing direction, and starts to align with the horizontal component of the magnetic field B after a few periods of walking.

In Fig. 3e, the member 10 first walks towards a glass capillary 14 mimicking a tunnel (inner diameter: 1 .62 mm). The tunnel impedes the walking gait of the member 10, such that it gets stuck at the tunnel entrance, because the cross-section of the glass capillary 14 does not allow the member 10 to lift up its front foot to proceed with further walking. To overcome this obstacle the member 10 then switches to the crawling mode of locomotion to pass through the tunnel. At this point the input is switched to a clockwise rotating magnetic field B, and the member 10 starts to crawl through the tunnel. After passing through the tunnel, the walking locomotion is resumed and the member 10 finally leaves the working space.

It has further been discovered that the magneto-elastic member 10 can additionally accomplish functional, shape change-enabled tasks. As an example, gripping and transportation of an object is shown in Fig. 4a. Fig. 4 shows various frames taken from a video depicting the different shapes the member 10 adopts while carrying out a pick up and place task. Fig. 4a shows how the member 10 approaches a cargo item 14 (nylon, 1 x0.8x1 .5 mm ³ ) by walking on a flat rigid surface, picks up the cargo by curling (frame V, |B| = 20 mT), transports away the cargo 14 by rolling and maintaining the curled shape, and releases the cargo 14 by uncurling at a new position.

Fig. 4b shows a dynamic and selective cargo 14 release. A paper tissue (0.5x0.5x0.1 mm ³ , used as an example of a model drug container) is bound to the member 10 by an extra appendage. After pre-bending the member 10 (frame I), B is quickly reversed (frame 11- IV) to open the appendage and release the cargo (frames V-VIII; cargo 14 indicated by the dashed line). The maximum magnetic field applied in these examples is |B| _max : 30 mT. The loading of the cargo 14 does not hinder the locomotive multimo- dality of the member 10. The scale bars shown relate to a scale of 1 mm. In order to produce a programmable magnetic soft composite member 10, a negative base mold is produced (not shown), for example of rectangular shape to form a rectangular shaped magnetized member. A decision is then made if the member 10 should have a varying magnitude of magnetization with varying direction of magnetiza- tion over its length I or whether it should have a constant magnitude of magnetization with varying direction of magnetization. If a constant magnitude of magnetization is desired then the base mold is filled with magnetizable material, i.e. an active component (actuatable component) is filled into the negative mold. An example of the active component can comprise a mixture of NdFeB and Ecoflex.

The member 10 comprising the active component is then installed on a jig and then subjected to an external magnetic field to create the different directions of magnetization of the segments s of the member 10. In order to permanently magnetize the member 10, the member 10 is subjected to a strong magnetic field having a strength of around 1 .5 T. This is typically achieved by placing and orienting the member and the jig within the magnetic field generated by the electric coils.

If a non-uniform magnitude of magnetization is required then the base mold is filled with a passive component in order to form a rectangular shaped member. For this pur- pose, a mixture of Al and Ecoflex is used as the passive component and simply poured into the mold in liquid form and allowed to cure.

Having calculated the ideal shape required to form a magnetized component of desired variable magnitude and orientation, a negative component mold for the magnet- ized component (also not shown) is formed in the passive component, for example by cutting out excess material of the passive component with a laser cutter to form a shaped band having a uniform or non-uniform width. The shape of the band of is based on a programmed magnetization profile. In a further step the magnetizable material, i.e. the active component, is filled into the negative mold formed by the passive component. An example of the active component can comprise a mixture of NdFeB and Ecoflex. The active component is poured into the component mold formed within the passive component and is cured in the component mold and thereby adheres to the passive component.

The member 10 comprising the passive component and the active component is then installed on a jig and subjected to an external magnetic field to create the different directions of magnetization of the segments s of the member 10. Again, the member 10 is subjected to a strong magnetic field having a strength of around 1 .5 T to permanent- ly magnetize the member 10. In particular the body of the magneto-elastic member 10 is fabricated by loading an Ecoflex-10 polymer matrix (density: 1 .04 g/cm ³ ) with neodymium-iron-boron micropar- ticles (MQP-15-7, Magnequench; average diameter: 5 μιτι, density: 7.61 g/cm ³ ) ac- cording to a mass ratio of 1 : 1 . The resulting soft magnetic elastomer has a density of 1 .9 g/cm ³ and a Young's modulus of 6.4x10 ⁴ Pa, as measured by a tensile testing machine (Instron 5943, Instron Inc.). During the fabrication, the pre-polymer is cast onto a flat poly(methyl methacrylate) plate to form a 240 μηη-thick film. After thermal polymer curing, the device is laser cut with the desired geometry.

Fig. 5a shows an example jig 16 that is formed by a glass rod. Other forms of jigs can be used to form specific profiles of magnetization. By using a non-uniform width negative component mold into which the active material is introduced a magnetic member 10 of non-uniform magnitude of magnetization is formed in which segments s of the member 10 also have varying directions of magnetization. In the present example, a magnetic field of 1 .5 T was applied to magnetize the member 10.

It should be noted in this connection that a magnetic field of at least 1 .2 T is required in order to saturate the magnetization of NdFeB.

A detailed description of the embedding process can be found in the publication by G. Z. Lum et al., titled "Shape-programmable Magnetic Soft Matter" and published in Proceedings of the National Academy of Sciences vol. 113, E6007-E6015 (2016), in particular in connection with Fig. 5 and the associated description of said publication. The publication by G. Z. Lum et al. also details the mathematical equations required to compute the shape of the negative component mold in order to form the member 10. European patent application EP16175341 .3 provides further details of this method of manufacture. The glass rod forming the jig 16 shown in Fig. 5a shows an elastomeric member 10 comprising NdFeB microparticles that is wrapped around the cylindrical glass tube (perimeter: 3.7 mm) and that is magnetized by a uniform 1 .5 T magnetic field produced by a vibrating sample magnetometer. The relative orientation _R defines the phase shift in the magnetization profile. Fig. 5b shows the programmed magnetization profile in the unwrapped member 10 of Fig. 5a, for _R = 45°.

The member 10 dimensions used in this work are shown in Fig. 5b, where I, w, and h are the length, width, and height of the member 10, respectively. To create the phase shift, β _κ , in the magnetization profile m (Fig. 5b), the rod is rotated by the angle β _κ during the member 10 magnetization. If not otherwise mentioned, a β _κ of 45° is assumed in the following. It should be noted that the native surface of the composite elastomeric material is hydrophobic and microscopically rough. Static advancing (1 16°±3°) and receding

(78°±2°) water contact angles were measured by the sessile droplet method through an automated goniometer routine that continuously increased and decreased the water droplet volume, respectively (Kruss DSA100). Surface roughness (R _a =0.63±0.02 μιτι and R _z = 4.37±0.28 μιτι) was measured by laser interferometry (Keyence VK-

X200). The members 10 surface was not subjected to any treatment before the experiments.

Fig. 6 shows an example of a six-coil C electromagnetic setup M to actuate the mag- neto-elastic member 10. The global coordinate system in X, Y and Z is indicated. The magnetic actuation setup is composed of three orthogonal pairs of custom-made electromagnets. The input currents driving the electromagnets are specified by software control signals through a custom electronic board. In the present setup a calibrated magnetic actuation matrix was used to map the control current input to the actual magnetic field in the working space. The effective working space along the X, Y and Z coordinates is restricted to 30 mm around the center of the reference system to exclude the effect of nonlinear components of the actuation matrix in close proximity of each coil. In particular Fig. 6 shows three pairs of coils C respectively arranged at a common axis, with the three common axes intersecting perpendicular to one another at the center of the apparatus. The spacing of each coil from the center is the same and in the specific example amounts to 37.5 mm. Through passing a current through the coils magnetic fields having a magnitude of 0 to 500 mT can be generated at the center of the arrangement. The magnetic field strength can be varied with a frequency of 1 kHz and the direction of the magnetic field can be varied with a frequency of 1 kHz.

In the global center of the working space, the relation between the current in coil and the magnetic field is:

in which are magnetic fields in +X, +Y,+Z directions and are in mT (millitesla); are the current flowing through the coil

pair, which contains two coils, in each direction and are in A (ampere). Positive current value is defined to produce positive magnetic field in this direction. For our system, the actuation matrix A is:

Based on the prescribed magnetization profile (Fig. 5b), a description will now be given to describe how the member 10 changes its shape when different external magnetic fields, B, are applied. Furthermore, the steering method is illustrated for the soft member 10. The analysis in this section provides a simple theoretical understanding for the experiments discussed in the following.

Quasi-static model

In order to discuss the quasi-static shapes that the soft member 10 can assume when subjected to B, the magnetization profile m, along the members 10 length, s, is assumed to be:

The variables m and <¾ represent the magnitude and spatial angular frequency of m(s), respectively, with Using the above m(s) and without any loss of gener-

ality, the members 10 local frame is used to describe the actuating magnetic field B (Fig. 5b). As the magnetic field B interacts with the magnetization profile, the soft member 10 experiences a locally varying stress that allows it to create desirable deflections. Fig. 7 shows a quasi-static analysis for the soft member 10. When the member 10 deforms, the bending moment acting on the infinitesimal element at steady state is Ads.

These static deflections can be described via the Euler-Bernoulli equation, where the rotational deflection along the member 10, 6(s), can be expressed explicitly by the actuating magnetic torque, _m , as:

where

The variables = hw, I = h ³ w /12 and E represent the member 10's cross-sectional area, second moment of area and Young's modulus, respectively, while R is the standard z-axis rotational matrix, given as:

Physically, the left side of Eq. S2 represents the magnetic actuation that deforms the member 10, while its right side represents the resultant deflections. Since magnetic forces are not applied onto the beam, the shear forces along the beam will be zero at steady-state and so can be neglected for the quasi-static analyses (see Fig. 7). For the locomotion modes that have small deflections, such as swimming and crawling, Eq. S2 can be simplified by approximating R with a simple 3 x 3 identity matrix, s≈ (along the body frame's x-axis), and where y refers to the vertical deflection. The gov-

erning equation then becomes:

Using Eq. S4, the shape of the member 10 can be described via the linear combination of two principal shapes, resembling respectively a cosine and a sine function. The first principal shape can be enforced when the direction of the magnetic field B is parallel to [cos β _κ sin /? _R 0] ^T , i.e., and as presented in (1 ), Eq. S4 can be double integrated to become:

In a similar manner, when B is parallel to [-sin(/? _R ) cos(/? _R ) 0] ^T , i.e.,

Eq. S4 can be double integrated to become:

While Eqs. S4-S6 can fully describe the achievable shapes for the swimming and crawling locomotion (see Fig. 1 b), they can only be used as an approximation for the other modes of locomotion, with the other modes of locomotion making use of large deflections. Despite the limitations, Eq. S5 can still be used to approximate large beam (member) deflections, since a large magnitude magnetic field B parallel to [cos /? _R sin /? _R 0] ^T creates a distorted cosine shape that resembles a semicircle (Fig. 1 b). While the cosine function cannot fully describe the semicircle-like shape, Eq. S5 suggests that a larger curvature can be achieved when the amplitude of the cosine function is increased. The capability to create the large curvature semicircle-like shape is essential to achieve most locomotion and transition modes.

Steering strategy

In the un-deformed state, the soft member 10 does not have a net magnetic moment. The member 10 can still produce a net magnetic moment upon deformation. As the net magnetic moment M _ne t of the member 10 tends to align with the applied magnetic field B, the direction of the magnetic field B can be controlled to steer the member 10 while it moves. M _ne t can be expressed

To achieve the multiple modes of locomotion, the magnetic field B is used to steer M _ne t and to allow the member 10 to move along an intended direction or to produce the necessary tilting (e.g., for walking).

In the following a description of the magnetic field B required to achieve each of the different modes of locomotion expressed by the shown member 10 will be given.

Moreover, additional experimental details necessary to combine all the modes of locomotion together will also be provided. Unless stated to the contrary the member 10 described in the following has a /3R of 45°.

Low-magnitude (< 5 mT) rotating B: Undulating swimming and Crawling

The undulating swimming and crawling modes of locomotion require similar undulating, time-varying shapes. Like the Taylor swimming sheet, the soft member 10 achieves its undulation by using a rotating magnetic field B of low magnitude and relatively high frequency.

Fig. 8 shows the rotating magnetic field B sequence enabling the undulating motion of the magneto-elastic member 10. Fig. 8a shows the required B plotted against time to undulate the member 10 to the right. The variables B _x and B _y represent the magnetic field along the x- and y-axis of the coordinate system, shown in Fig. 8b. Fig. 8b shows the two principle shapes as described by Eq. S5-S6 (reproduced from frame I in Fig. 1 b). The scale bars shown scale to 1 mm. Based on Eq. S5-S6, when the field angle of the magnetic field B, represented by a, is 135° or 225° (Fig. 8b), a small magnetic field |B| (3 mT) enables the member 10 to produce a shape that resembles a cosine or a sine function, respectively. As a is subsequently rotated, the member 10 produces a traveling wave along its body. The deformation produced by the time-dependent undulation can be described mathematically as:

where the variable T represents the period of the rotating magnetic field B.

High magnitude (> 7 mT), low frequency (< 15 Hz) B: Rolling, Immersion, Meniscus Climbing and Landing

Fig. 9 shows the rotating sequence of the magnetic field B enabling the rolling mode of locomotion. Fig. 9a shows the required B plotted against time (normalized to one period) to roll clockwise. The variables B _x and B _y represent the magnetic field along the x- and y-axis of the coordinate system, shown in Fig. 9b. Fig. 9b shows the semi-circle shape required for the rolling locomotion (reproduced from frame I in Fig. 2e). The scale bars indicated an actual size of 1 mm.

To perform the rolling locomotion shown e.g. in Fig. 2e, a large magnitude magnetic field B (15 mT) is first applied along a field angle of 157° to deflect the member 10 into a distorted cosine shape resembling a semicircle-like shape (see Fig. 9b). The deflected member 10 produces a M _ne t described by Eq. S7, allowing the member 10 to be steered via a gradual rotation of the direction of the magnetic field B. During this period of time, a rigid-body magnetic torque, Xrigid-body, is induced by the member 10, which can be described mathematically as:

Once the rotation of the magnetic field B is stopped, the M _ne t Of the member 10 will be aligned with its final direction, and T _rigid _ _body will be reduced to a null vector.

Fig. 10 shows the sequence of the magnetic field B to perform the meniscus climbing mode of locomotion. Fig. 10a shows the required magnetic field B plotted against time to climb the meniscus on the right. The variables B _x and B _y represent the magnetic field along the x- and y-axis of the coordinate system, shown in Fig. 10b. When the magnetic field B has a field angle of 315°, the member 10 produces a shape that resembles a semi-circle. The member 10 pins on the water surface with its feet under deformation. The scale bar shown depicts an actual size of 1 mm. This rigid-body torque is necessary for meniscus climbing. The required magnetic field B sequences to perform meniscus climbing are shown in Fig. 10a. The initial magnetic field B (|B|=1 1 .7 mT, a = 315°) is imposed to let the member 10 assume a vertically symmetric pose with upward curvature (Fig. 10b). In 400 ms the member 10 achieves a stable position along the meniscus known to depend on both the initial meniscus profile and the curvature of the device. The magnetic field B is subsequently rotated clockwise toward a = 157° to progressively displace the position of the member 10 along the meniscus and raise its center of mass. Once the device reaches contact with the solid substrate, the magnetic field B is turned off and the member 10 rests on the edge of the adjacent platform (see also Frame I in Fig. 2d).

Fig. 1 1 shows the rotating B sequence required to perform immersion. Fig. 1 1 a shows the required B plotted against time (normalized to one period) to rotate counterclockwise and immerse into the water. The variables B _x and B _y represent the magnetic field along the x- and y-axis of the coordinate system, shown in Fig. 1 1 b. Fig. 1 1 b shows the initial curving of the member 10 for immersion (reproduced from frame I in Fig. 2c). The scale bar shown depicts an actual size of 1 mm.

The sequence of the magnetic field B for immersion is shown in Fig. 1 1 . In this mode of locomotion, the member 10 rotates counterclockwise to reduce the contact with the water-air interface to a minimum, right before a quick flip in the direction of B intervenes to disengage it and let it sink (0.9 time unit in Fig. 1 1 a).

Fig. 12 shows the rotating B sequence to perform landing. Fig. 12a shows the required magnetic field B plotted against time to escape the capillary attraction of the water. The variables B _x and B _y represent the magnetic field along the x- and y-axis of the coordinate system, shown in 12b. Fig. 12b shows the initial curving of the member 10 for landing (reproduced from the frame I of Fig. 2d). During the landing mode of locomotion water can easily dewet the surface of the member 10 because its high receding contact angle leaves the receding contact lines unpinned. The magnitude of the applied magnetic field B decreases at the 0.2 time unit in Fig. 12a from 1 1 .7 mT to 5 mT. Before this time point, most of the member 10 has peeled away from the water surface and a large magnetic field B is no longer nec- essary.

Fig. 13 shows the rocking sequence of the magnetic field B required to perform the walking mode of locomotion. Fig. 13a shows the required magnetic field B plotted against time to walk rightward. The actuation is divided into four phases, correspond- ing to Frame I -V in Fig. 2f. Note that the actuation is periodic, so that Frame I corre- sponds to Frame V. Fig. 13b shows how when the magnetic field B has a field angle of 315° the member 10 produces a shape that resembles a semi-circle (reproduced from the frame I of Fig. 2f). The scale bar shown depicts an actual size of 1 mm. Fig. 14 shows the flipping magnetic field B sequence enabling the jelly fish-like swimming motion mode of locomotion of the magneto-elastic member 10. Fig. 14a shows the required B plotted against time to make the member 10 swim like a jellyfish upwards inside water. The variable B _y represents the magnetic field along y-axis of the coordinate system, shown in Fig. 14b. B _x is zero throughout the period. Fig. 14b shows the two shapes assumed by the member 10 during swimming. The scale bar shown depicts an actual size of 1 mm.

The jellyfish-like swimming shown in Fig. 14 in bulk liquid is implemented by periodic flipping of the magnetic field B along the vertical direction to give a fast power (down- ward) stroke and slow recovery stroke. In this way the member 10, functioning at a Reynolds number Re -2.3, is shown to swim upward against the gravitational field, and even against the downward magnetic field gradient eventually present at the vicinity of the bottom coil of the working space. Flipping the magnetic field B can also make the member 10 jump. To this end Fig. 15 shows snapshots of frames taken to illustrate straight and directional jumping by way of the magneto-elastic member 10. The first six video frames and time-lapse (last column) are shown for each jumping sequence. The same member 10 with _R = -90° was used in all cases. The recording speed was adapted to the dynamics of each jump type. Fig. 15a shows Type-1 straight jumping (550 fps). There is almost no net magnetization in the second frame as the member 10 is almost flat. In this example the maximum magnetic field B _max is12 mT.

Fig. 15b shows Type-2 straight jumping (355 fps). In this example the maximum mag- netic field B _max is 18.9 mT. Fig. 15c shows directional jumping (720 fps), which is also shown in Fig. 2e. The direction of the magnetic field B (18.9 mT) applied is switched from 189° to 0° and maintained for 25 ms. In the frames where the magnetic field B is applied, the instantaneous net magnetization vector of the member 10 is indicated by the arrow. Fig. 15d shows the magnetization profile used for the member 10 in (a-c), and its magnetic response under large |B|. The scale bar shown depicts an actual size of 1 mm.

All types of jumping are achieved by quickly flipping the direction of the magnetic field B via a step or a fast ramp function. In order to effect this a single pair of coils is pref- erably used to create a (vertical) magnetic field B along the Y-axis of the global reference system and the factors affecting the jumping performance.

When a R = 45° is used, a magnetic field B along the (1 ,1 ) direction in the X-Y plane can be used, though this is prone to more errors since two pairs of coils have to be used.

In Type-1 straight jumping, the member 10 is first pre-bent by the magnetic field to curl upward and stand on its center (frame I of Fig. 15a). The direction of the magnetic field B is then quickly reversed to force the member 10 to flatten out (frames II and III of Fig. 15a) and hit the substrate with its feet. An example of the modeling framework used to optimize this type of jump is shown in the following. In contrast to this, Type-2 straight jumping does not require pre-bending, as the magnetic field B is only applied to curl the member 10 to jump (Fig. 15b). This difference between Type-2 and Type-1 jumping is characterized by using: 1 ) different ramping time of external magnetic field B; 2) different magnitude of the magnetic field B used; 3) a different magnetic field B spatial gradient; and 4) for two types of m.

Directional jumping (Fig. 15c) is also featured in the multimodal locomotion experiment shown in Fig. 3b. The results shown in Fig. 3e are re-plotted in Fig. 15c with more details. For this type of jump, the member 10 is initially deformed into a shape that resembles a cosine function. Subsequently, the direction of the applied magnetic field B is slightly rotated (-9°) about the longitudinal axis of the member 10 to specify the rotational direction during the jump.

Following this a step change in the direction of the magnetic field B is imparted, this simultaneously incepts both rigid-body rotation (clockwise here, to align the net magnetization to the external magnetic field B direction), and shape change (frames l-lll in Fig. 15c). During the rotation, the shape of the member 10 changes to an inverted co- sine, increasing the momentum of the member 10 before it strikes the substrate (frames lll-V in Fig. 15c).

It should be noted in this connection that the jumping mode of locomotion of the member 10 is sensitive to the initial pose and shape of the member 10. Exact consecutive reproduction of the same jump may be hindered by the inherent viscoelasticity of the body of the member 10 and by the associated shape hysteresis.

In the following a discussion will also be presented on how the dimensions of the member 10 affect the performance of each mode of locomotion, in particular how the dimensions of the member 10 affect its swimming and crawling speed. In order to facil- itate this discussion, simple theoretical models are employed that can approximate the physics of each mode of locomotion. The derivation of these models is a first step towards understanding the full physics of a miniature soft member 10 with multimodal modes of locomotion.

The speed of the swimming and crawling modes of locomotion depends on the highest rotational frequency of the magnetic field B that the member 10 can respond to, the mechanical bandwidth of the member 10, i.e. its non-zero fundamental natural frequency can be computed. This computation can be achieved by using a simple vibra- tional analysis where small deflections for the member 10 are assumed. Using this assumption, Newton's second law is initially used to analyze an arbitrary infinitesimal element of the member 10, dx, at a time t. Fig. 16 shows a vibrational analysis of the member 10 by showing the forces and torques acting on an infinitesimal element of the member 10.

Based on Newton's force law, the equation of motion of such an element along the vertical axis is given as:

where C, v and p represent the damping coefficient, shear force and density along the member 10's body, respectively. In Eq. S10, -Cwy and pAy represent the damping and inertia forces acting on the element, respectively.

In a similar manner, τ = Ja can be used to obtain the rotational equation of motion for this element about its bending axis:

The variable M _b represents the bending moment acting on the element, and based on the Euler-Bernoulli equation, it is also equivalent to Therefore, by substituting

Eq. S1 1 into Eq. S10, one can obtain the governing equation for the beam:

The well-established mode-shape analysis is used to evaluate this partial differential equation. A free-vibration analysis is initially used to obtain the mode shapes of the member 10, i.e. the damping effects and magnetic actuation are temporarily removed from Eq. S12. Subsequently, the variable separation method for y is used such that

where F(x) is solely a function of x and G(t) is only a function of time. By substituting Eq. S13 into Eq. S12 (without the damping effects and magnetic actuation) and rearranging the terms, one obtains:

Based on classical vibration analysis, both sides of Eq. S14 are independent of x and t, and they can be equated to become a function of the R ^th natural frequency, of the member 10 (4):

After solving the homogeneous equation in Eq. S16, one can obtain both ω _{η Κ} ² and the R ^th mode shape of the member 10 (F _R ):

where

and the variables A _R , B _R , CR, and D _R are constants dependent on the boundary conditions of the member 10. Because the member 10 has free ends, the following free- free boundary conditions are substituted into Eq. S17:

Physically, the boundary conditions in Eq. S18 dictate that there are no shear forces and bending moments acting on the free ends of the member 10. Based on these boundary conditions, one can obtain the following characteristic equation for the member 10:

Once the characteristic equation has been solved numerically, the two lowest natural frequencies of the member 10 are given as:

The first mode, ω _{η 1} , represents the rigid-body motion of the member 10, while the second mode shows the non-zero fundamental natural frequency of the member 10. As the first mode does not describe the shape change of the member 10 with respect to the actuating magnetic field's frequency, the presented frequency response analysis begins with the second lowest natural frequency, i.e.s ω _{η 2} . Based on Eq. S20, this natural frequency can be expressed as:

Using Eq. S21 and the Taylor swimming sheet model, one can approximate the achievable swimming speed, v _SW im, of the member 10. Although the selected fluidic boundary conditions result in a slight deviation from the Taylor model, i.e. the presented member 10 swims on the water surface and not far away from walls, it is assumed that the following equation can still capture the fundamental dynamics of the swimming locomotion:

where f represents the frequency of the traveling wave (coincident with the rotation frequency of B). As v _SW im can be maximized by applying the highest f and since the highest f is limited by one can rewrite Eq. S20 as:

Based on Eq. S23, the swimming speed of the member 10 can be increased when its I and h are increased and decreased, respectively. Similar to the swimming mode of locomotion, the member 10 can crawl faster when it can produce a faster undulating frequency. However, in order to enhance the crawling speed one would have to increase the second lowest natural frequency of the member 10, and Eq. S21 implies, in contrast to Eq. S23, that /i should be increased while / should be decreased.

The rolling mode of locomotion also uses a rotating magnetic field B to control the soft member 10. However, in this case the magnetic field B rotation frequency f has to be much slower (/ < 15 Hz in the present example) than the one utilized in the swimming and crawling mode of locomotion to prevent the member 10 from creating a traveling wave along its body. Therefore, it will not be necessary to use vibrational analyses to evaluate the performance of this mode of locomotion.

It is assumed that the rolling mode of locomotion is similar to that of a wheel, such that it is more efficient when the member 10 becomes rounder. A rounder shape will re- quire the soft member 10 to produce deformations that have large curvatures. While the cosine function in Eq. S5 cannot fully describe the semicircular shape of the mem- ber 10 for the rolling mode of locomotion, this equation suggests that a larger curva- ture can be achieved when the amplitude of this cosine function, i.e., can be

increased. As m, B, E and other constants are independent of the dimensions of the member 10, one can conclude that a larger curvature can be easily accomplished when is maximized. Since larger curvatures can be achieved by

increasing / and decreasing h.

Similar to the rolling mode of locomotion, the effectiveness of meniscus climbing is highly dependent on the maximum achievable curvature of the member 10. As a re- suit, one can use Eq. S5 and the above discussion on rolling to approximate how the dimensions of the member 10 can affect the climbing behavior of the member 10. This equation suggests that a larger curvature can be achieved by increasing / and decreasing h. Jumping

Fig. 17 shows the boundary condition for the Type-1 straight jump (see also Fig. 15a). The analysis can be simplified by considering a fixed-free half beam, since the center point of the member 10 neither rotates nor translates and the motion is symmetrical about this point.

As it is very challenging to derive a mathematical model that can accurately describe how the dimensions of the member 10 can affect its jumping performance, only a simple vibrational analysis will be provided as a first step towards understanding this mode of locomotion. The presented analysis is performed to approximate the simplest jumping mode - Type-1 straight jump in Fig. 15a. To proceed with this vibrational analysis, one assumes that the shape of the member 10 begins with a cosine function in stage I (Fig. 17a). Furthermore, after the member 10 is subjected to a step change in B direction (i.e. flipping B for 180°), the shape of the member 10 will tend to change into an inverted cosine function. During this transition, the member 10 experiences an impact with the substrate, allowing itself to jump (stage II).

The formulation can be reduced to analyze only half of the beam because the jumping motion is symmetrical about the center point of the member 10, which neither rotates nor translates, and its boundary condition can be simplified into a fixed-free type (Fig. 17).

In order to facilitate the formulation in the following, one uses V to represent 1/2 and x' = x - V. Using these variables, one can mathematically represent the boundary conditions for half of the member 10 as:

Although the boundary conditions for the jumping mode of locomotion is different from the swimming mode of locomotion, the vibrational analysis from Eq. S10-17 is still valid. By substituting the boundary conditions in Eq. S24 into Eq. S17, one obtains the following characteristic equation:

Once the R ^th natural frequency in Eq. S25 has been solved numerically, this value is substituted back into Eq. S17 to obtain the corresponding mode shape:

F _R (x)

By substituting the obtained mode shapes into Eq. S12, rearranging the terms and normalizing the inertia terms, the following equation is obtained:

Because the magnetic actuation on the right side of Eq. S27 can be computed as:

To simplify the discussion, only the fundamental natural frequency and mode shape are used to describe the jumping mode of locomotion. While the exclusion of other mode shapes might introduce errors into this prediction, it is believed that the discus- sion is still valid, since the fundamental mode shape is the most critical component in Eq. S28. Therefore, the orthogonal property for vibration is used and Eq. S28 is multiplied with F ₁ and subsequently both sides are integrated over x' from x' = 0 to ' = V to obtain the following equation:

The characteristic equation in Eq. S25 numerically evaluates that and

by using this value one numerically determines

. By substituting these values back into

Eq. S27 and normalizing the inertia term, one obtains:

At stage I, the member 10's initial shape can be expressed as

where A ₀ is a positive real number that represents the amplitude of the cosine function. Similar to Eq. S27, one can use the orthogonal property of classical vibration analysis to extract the initial condition for G^t = 0), i.e., one can multiply the initial shape of the member 10 with _x and subsequently integrate the new function over x ' from 0 to I ' . Based on numerical evaluation, one obtains G (t = 0) = G _{1 Q} l, where Gi _, o is a positive real number. Because the member 10 does not have an initial speed, one can conclude that the second initial condition can be expressed as G^t = 0) = 0. Therefore, by using these two initial boundary conditions and assuming the system is under- damped (i.e. low magnitude of C) one can solve Eq. S30 to obtain:

Where

The kinetic energy for the member 10 before it impacts the substrate at time i _f , E _kEi£ , can then be approximated numerically as: where V represents the volume of the member 10. After striking the substrate, Newton's law of restitution dictates that the kinetic energy of the member 10, will be

come:

where / represents the coefficient of restitution and it is a positive real number that is between 0 and 1 . One assumes that can be converted fully to gravitational poten

tial energy such that:

where z is the highest point of the jump. Therefore, by substituting Eq. S32 and S33 into S34, and rearranging the equation, one can approximate z as:

To maximize z, the magnitude of and y ² have to be maximized. However,

as is highly dependent on other parameters like C,p, E and if, it is difficult

to deduce the effects of h and I for such a complex function. Likewise, it is also difficult to deduce the mathematical relationship between y ² and the dimensions of the member 10, i.e., its w, h and /. As a result, it is believed to be necessary to perform a numerical optimization to determine the optimal values o h, w and I such that the mem- ber 10 can maximize z.

Walking and Jellyfish-like Swimming

The net stride of the member 10 in one walking cycle is dependent on the difference between the curvature of the member 10 in Steps 3 and 4. Thus, when a larger curvature can be achieved in Step 4, the member 10 can produce a bigger stride. By using similar analyses as presented in the foregoing, one can deduce that a larger stride can be achieved per cycle when / is increased and h is decreased. On the other hand, the walking speed also depends on the bandwidth of the member 10, i.e. the capability to respond to high frequency periodic magnetic fields B. In this instance the bandwidth of the member 10 is approximated with its fundamental natural frequency, so that based on Eq. S21 the bandwidth can be increased when / is decreased and h is increased. Because the criteria to produce a larger stride are contrary to the criteria for higher bandwidth, it is possible to perform numerical optimization methods to determine the optimal dimensions of the member 10.

In a similar manner, greater swimming propulsion can be generated in one swimming cycle (jellyfish-like swimming) when the speed difference between the power and re- covery strokes is increased. The fastest achievable speed is dictated by the bandwidth of the member 10. Therefore, it is beneficial to decrease I and increase h such that the fundamental natural frequency of the member 10 can be maximized.

However, decreasing / and increasing h will increase the overall stiffness of the mem- ber 10, making it harder to create a large stroke length that is necessary for generating large propulsions. Because the criteria to produce a larger stroke is contrary to the criteria for higher bandwidth, it is assumed that a numerical optimization is necessary to determine the optimal dimensions of the member 10. Table 2 summarizes the design considerations discussed above for the member 10 under the control of a magnetic coil setup with negligible settling time.

Given the magnetic actuation, it is possible to use the spatial gradients of the magnetic field B to exert magnetic pulling on the member 10 to enhance its locomotion. For in- stance, after the member 10 has jumped, the deformed member 10 has a net magnetic moment that allows the use of gradient pulling to modulate jump height and distance.

Theoretically, magnetic gradient pulling may also levitate a magnetic member 10 and make it traverse across different terrains and obstacles. However, this approach is not practical. First, in absence of closed loop feedback control, the trajectories generated by gradient pulling are generally jerkier than those achievable by the proposed locomotion. Second, the required spatial gradients of the magnetic field B for the member 10 to levitate across all obstacles and terrains, estimated to be 3.8 T/m, are difficult to obtain. Moreover, when envisioning the use of the member 10 in healthcare, high spatial gradients of the magnetic field B are considered adverse for e.g. medical application, because they can create undesirable forces on implanted objects within a patient's body, e.g. an implantable cardioverter defibrillator (ICD). Although advancements in the design of ICDs are underway, the shock coils of many ICDs are ferro- magnetic and they can experience significant adverse translational motions when they are subjected to such high spatial gradients. In contrast, although the ICD devices will also be affected by the required magnetic actuation, most ICD devices can fully function when the magnitude of the magnetic field B is lower than 500 mT, and this value is well above the magnetic field strength utilized for the soft member 10 (-20 mT) pre- sented herein.

A design of a modified member 10 for selective, magnetically triggered drug release is shown in Fig. 18. To this end Fig. 18a shows that an extra strap 18 is added to the body of the original member 10. For clarity, only the main component of magnetization on the extra strap 18 is marked out in Fig. 18a, as it is the dominant and functional component in releasing a drug that is provided as the cargo 14. The magnetization of the member 10 main body is the same as that shown in Fig. 5b.

During locomotion, the cargo 14 is mechanically bound on the member 10 by inserting the strap head 20 into a hole 22 on the body 24 (Fig. 18b). When a large magnetic field B is applied along the y-axis, B _y (-19 mT), then the magnetic torque bends and unlocks the strap to release the cargo (Fig. 18c). Note that other components of magnetic field (Βχ, B _z ) may also be present. Due to the restriction of the hole 22, a smaller By (<18 mT in our case) is not sufficient to open the strap 18, which makes normal lo- comotion possible. The design may be optimized to also carry liquid or powder based loads as cargo 14. in the foregoing approximate quasi-static and dynamic models have been shown to describe how the width w, length /, and thickness h of the member 10 can be designed to optimize each mode of locomotion (see also Table 2). Hence observations were made that w only affects jumping and jellyfish-like swimming, whereas / and h are subject to conflicting requirements for different locomotion modes. Improved meniscus climbing and rolling is predicted when / is increased and h is reduced, as this increas- es the member's bending compliance Conversely, fast and efficient crawling

requires limited member compliance to preserve its ability to track high frequency control signals. The walking gait imposes by itself conflicting requirements on / and h. A more compliant member body is preferred to produce a larger net stride and larger body curvatures, which implies the increase and decrease of / and h, respectively. However, this would also effectively decrease the bandwidth of the member, reducing the maximum frequency of steps in the walking gait. The present study only elicits a feasible design that enables all reported modes of locomotion.

Table 2 shows design considerations of the member 10 for multimodal locomotion. The arrows respectively indicate whether an increase or decrease in the length, height or width have beneficial effects on which type of locomotion in comparison to the member 10 discussed in the foregoing, so as to optimize it for a certain mode of locomotion.

The magneto-elastic, rectangular sheet-shaped, soft robot 10 is made of silicone elastomer (Ecoflex 00-10) embedded with hard magnetic neodymium-iron-boron (NdFeB) microparticles that have an average diameter of 5 μηη. The surfaces of the robot are hydrophobic, and they can potentially be made biocompatible (SI section S1 C). By following the magnetization process described in SI section S1A, the robot can be programmed to have a single-wavelength harmonic magnetization profile m along its body (Fig. 19a and Fig. S1 ). After m is programmed, the robot can be controlled by a time-varying magnetic field, B, to generate different modes of locomotion. Unless oth- erwise specified, B is spatially uniform, and therefore no magnetic forces are applied to translate the robot (SI section S15). The uniform B, however, can control the robot's morphology and steer it to move along a desired direction. To describe the effects of

B, we express with respect to the robot's body frame (Fig. 19a) where

B _xy represents the xy plane components of B, i.e., The interaction between B _xy and m produces spatially-varying magnetic torques that deform the robot, and hence controlling B _xy allows us to generate the desired time-varying shapes for the robot. As the deformed robot possesses an effective magnetic moment M _ne t (Fig. 19b), which tends to align with B, we can control B _z to rotate the robot about its y-axis, steering it along a desired direction (see SI section S3B(II)).

Depending on the magnitude of B _xy , i.e., B _xy , the robot exhibits different shape- changing mechanisms (Fig. 19b and SI section S3A-B). When B _xy is small {e.g., B _xy ≤ 5 niT) and B _xy is aligned along the two principal directions shown in Fig. 19b(ll-lll), the prescribed m produces a sine or a cosine shape for the robot. Because the robot's deformation is small in such conditions, orienting B _xy away from the principal directions generates a weighted superposition of the two basic configurations. Thus, we can create a traveling wave along the robot's body by using a rotating B _xy that has a small constant magnitude. As the robot's M _ne t is always parallel to the applied B _xy in small- deflection conditions, the robot does not experience any rigid-body magnetic torque and consequent rotation about its z-axis (SI section S3B(I)). Conversely, when B _xy has high magnitude {e.g., B _xy = 20 mT) and is aligned along the principal axis shown in Fig. 19b(IV-V), the robot undergoes a large deflection shape-change to deform into either a 'C- or a 'V'-shape. However, if the direction of B _xy is not along this principal axis, the deformed robot generates a large M _ne t that is generally non-parallel to the applied field, and this makes the robot rotate about its z-axis until its M _ne t aligns with B _xy (Fig. 19c and SI section S3B(I)). At the end of this rotation, the robot will assume its 'C- or 'V'-shape configuration because the generated M _ne t in these configurations is naturally aligned with the applied B _xy . Using this mechanism, we can control the robot's angular displacement about its z-axis to enable locomotion modalities like rolling, walking and jumping.

By using the steering and shape-changing mechanisms, we demonstrate all of our robot's locomotion modes in Figs. 2 and 3. When completely immersed in water, the robot can swim upwards and overcome gravity (Fig. 20a, Supporting and SI section S10). A periodic B with time-varying magnitude along the principal axis allows the shape of the robot to alternate between the 'C- and 'V'-shapes, enacting a gait similar to jellyfish swimming. Inertial effects at Reynolds number ranging from 74 to 190 permit this time-symmetrical but speed-asymmetrical swimming gait to produce fluid vorti- ces that propel the robot to reach the water surface (Fig. 3c, and SI Fig. S37). Upon emersion, the soft robot strongly pins at the water-air interface by exposing its hydrophobic surface to air.

Inspired by beetle larva that overcome frictionless barriers by performing quasi- static work on liquid-air interfaces, the soft robot can climb up a water meniscus by deforming into a 'C'-shape to enhance its liquid buoyancy without extra energy expenditure (Fig. 20b, and SI section S7). Upon meniscus climbing and reaching contact with an adjacent solid platform, a slow rotating B will make the 'C'-shape robot rotate about its z-axis. The hydrophobicity of its surface allows the robot to be peeled away from the water surface by such rotation (Fig. 20c, and SI section S1 1 B). In contrast to meniscus climbing, the robot can also dive into the liquid bulk by disengaging from the water-air interface via a fast sequence of downward bending, rotation and flipping (Fig. 20d and SI section S1 1A).

In nature, soft-bodied caterpillars use rolling locomotion to escape from their predators, because this is an efficient and fast way to sweep across solid terrains. In- spired by caterpillars, our robots can also roll directionally over a rigid substrate or dive from a solid onto a liquid surface (Figs. 20e and 3a). This locomotion is enabled by a high-magnitude rotating B {e.g., B = 18.5 mT), which allows the robot to roll in its 'C'- shape configuration (SI section S5). However, the curled-up robot cannot roll across substrate gaps wider than its diameter, which can instead be accomplished by walk- ing.

Walking is a particularly robust way to move over unstructured surfaces and affords precise tuning of stride length and frequency (Fig. 20f, and SI section S6). Inspired by the inchworms' walking gait, the robot can walk along a desired direction when we use a periodic B to sequentially adapt its tilting angle and curvature. In each walking cycle, the robot first anchors on its front end to tilt forward so that it can pull its back end forward. The robot then anchors on its back end to tilt backwards and extends its front end to achieve a positive stride in a single cycle.

When the walking robot is blocked by narrow openings, it can mimic another caterpillar locomotion and use an undulating gait to crawl through the obstacle (Fig. 20g, and SI section S9). Crawling is encoded by a rotating B to produce a longitudinal traveling wave that propels the robot along the wave's direction. A similar control sequence additionally enables the robot to produce an undulating gait to swim efficiently on liquid surfaces (Fig. 3a). In contrast to crawling, however, the undulating swimming direction is antiparallel with the traveling waves' direction. While previous robots with multi-wavelength, harmonic magnetization profiles can also enable the undulating swimming locomotion ⁷ , such robots have not been able to create the critical 'C and 'V'-shapes necessary to realize multimodal locomotion.

Similar to nematodes, the soft robot can jump over obstacles, which are too high or time consuming to roll or walk over, by imparting an impulsive impact on a rigid surface (Fig. 20h, and SI section S4). The B control sequence prompts both the robot's rigid-body rotation, which specifies the jumping direction, and elastic deformations to maximize the momentum of its free ends before striking the substrate. This sequence of B is specified in the robot's local yz plane where B _y is used for inducing the shape-changing mechanism, while the rigid-body rotation of the robot is induced by both B _y and B _z .

To illustrate the robot's potential to navigate across unstructured environments (SI section S13), we demonstrate that the robot can use a series of locomotion modes to fully explore a hybrid liquid-solid environment (Fig. 3) and a surgical human stomach phantom (Fig. 21 a, and SI section S14A). Heading toward an in vivo ultrasound- guided operation, we also show that the robot can be visualized by an ultrasound medical imaging device as it rolls within the concealed areas of an ex vivo chicken muscle tissue (Fig. 21 b, SI section S14B, and Fig. S44). The soft robot can additionally accomplish functional tasks like gripping an object and transporting it to a targeted location (Fig. 21 c), as well as ejecting a cargo that is strapped on the robot (Fig. 21 dand SI section S14C).

In addition to magnetic field-induced torques, magnetic gradient-based pulling forces could also be used to enhance locomotion performance (e.g., speed, jumping height). Moving along this direction, we show that the jumping height can be increased by adding magnetic gradient-based pulling forces, and we will explore other similar possibilities in the future. Using gradient-based pulling exclusively may however be detrimental, as the dynamics of this actuation method is inherently unstable. From a practical standpoint, gradient-based pulling methods are also less energy efficient than locomotion propelled via magnetic field-based torques (SI section S15).

The lack of an on-board actuation method prevents the proposed robot from operating in large open spaces, making it unsuitable for outdoor applications like envi- ronment exploration and monitoring. Furthermore, the current demolding process creates a pre-stress in the magneto-elastic material that induces a small residual curvature in the robot when it is in rest state (Fig. 19b(l) and SI section S1 D). While the pre- stress does not hinder the robot from achieving multiple modes of locomotion and could be reduced through improved fabrication, it induces small errors in the predicted robot shapes (Fig. 19b(ll, V)) and partly affects the model matching of the experimental data for the walking and undulating swimming speeds (SI sections S6 and S8).

To better understand small-scale soft-bodied locomotion, we have devised theoretical models to perform a scaling analysis on how the robot's dimensions (L, w, and h, shown in Fig. 19a) would affect the jumping, rolling, walking, meniscus-climbing and undulating swimming locomotion modalities (SI sections S4-8). As the theoretical models for the crawling and jellyfish-like swimming locomotion are too difficult to be derived, we have instead used the experimental data in SI sections S9-10 to derive corresponding fitting models. Based on our theoretical and fitting models, it is predict- ed that a larger L and a smaller h are always preferred for multimodal locomotion as it can help the robot to move faster and jump higher. The models also suggest that w would only affect the jellyfish-like swimming locomotion and minimizing w would increase the swimming speed. There are, however, practical lower bounds for both h and w, because our current fabrication technique has difficulties creating robots that have h < 40 μηη and w < 0.3 mm. Likewise, the practical upper bound of L is typically constrained by the size requirements of specific applications and the maximum allowable workspace of the electromagnetic coil setup that generates the spatially uniform B. A more detailed summary for the scaling analysis and fabrication limits can be found in SI section S12.

To validate the theoretical models, we compared them against extensive experimental characterizations conducted across robots with differing dimensions. In general, except for the undulating swimming locomotion, the experimental data agree well with our models (see SI sections S4-8, SI section S12, and Table S4). Detailed discussions pertaining to the theoretical and experimental discrepancy for the undulating swimming locomotion can be found in SI section S8. While the present study only elicits a feasible robot design that enables all reported modes of locomotion, these analyses can also provide useful design guidelines for optimizing the performance of future miniature robots that have multimodal locomotion.

In the future, we will also use our proposed robot to study small-scale soft- bodied locomotion on other complex terrains such as within non-Newtonian fluids or on granular media. Moreover, we will scale down the proposed robots to sub- millimeter size scale and demonstrate their potential in vivo medical applications. Data availability statement

All data generated or analyzed during this study are included in the published article and its supplementary information files. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reason- able request.

S1 - Robot fabrication and characterization

In this section, we describe the robot's fabrication process and mechanical properties (section S1A), its surface properties (section S1 B), how it can potentially be made biocompatible (section S1 C), and also its residual strain energy (section S1 D).

A. Fabrication process and mechanical properties of the robot

The base material of the magneto-elastic robot body is an Ecoflex 00-10 polymer matrix (Smooth-On Inc.; density: 1 .04 g/cm ³ ) loaded with neodymium-iron-boron (NdFeB) magnetic microparticles (MQP-15-7, Magnequench; average diameter: 5 μιτι, density: 7.61 g/cm ³ ) according to a mass ratio of 1 :1 (the ratio of the mass of NdFeB microparticles to the mass of Ecoflex-10 elastomer). The NdFeB microparticles' high remanent magnetization enables us to generate effective magnetic actuation for the robots while their high magnetic coercivity ensures that our robot's magnetization profile is maintained during magnetic actuation: NdFeB particles requires more than 600 mT to be demagnetized, and our maximum B is limited to 50 mT. The resulting magnetic elastomer has a density (p _r ) of 1 .86 g/cm ³ _.

To fabricate the robots, the pre-polymer is cast onto a flat poly(methyl methacry- late) plate to form a 185 μηη-thick film. After thermal polymer curing, the device is demolded by using the laser to cut it into the desired geometry before its removal from the substrate plate. The robot dimensions are shown in Fig. S1 b, where L, w and h represent the length, width, and height of the robot, respectively. Unless specified otherwise, L = 3.7 mm, w = 1 .5 mm, and h = 185 μιτι are assumed. These are the dimensions for the robot described e.g. in Fig. 19. The harmonic magnetization profile of the robot is implemented by wrapping its soft body around a cylindrical glass rod with a circumference of 3.7 mm and subjecting it to a large, uniform magnetizing field of 1 .65 T (Fig. S1 a). To create the phase shift, β _κ , in the magnetization profile m (Fig. S1 b), the rod is maintained by an angle of /? _R during the robot's magnetization process. Unless specified otherwise, β _κ is assumed to be 45°. We experimentally evaluated the magnitude of m to be 62,000 ± 10000 A/m for robots that have a 1 :1 mass ratio (Fig. S1 c).

The effective Young's modulus (E) of robots, which has a 1 :1 mass ratio, is experimentally evaluated to be 8.45x10 ⁴ ± 2.5x10 ³ Pa through a tabletop test system (MTS Nano Biomix, Agilent Systems) using the accompanying software (Nano Suite V5, Agilent Systems). Figure S1 c-d additionally summarizes the characterization of |m| and E, respectively, for similar magnetic elastomers with different mass ratios. The mean and standard deviation of each data point in these sub-figures are computed from the experimental results of three samples. According to these experimental data, increasing the content of magnetic particles in the polymer matrix increases |m| more signifi- cantly than it increases E. This suggests that robots with higher NdFeB concentration (larger mass ratios) can deform more easily, as they can generate higher magnetic torques, while becoming only slightly stiffer. From the fabrication perspective, it is difficult to further increase the NdFeB mass content beyond the 1 :1 mass ratio as the base material Ecoflex-10 no longer cures well. Hence, all the robots used in this work are fabricated based on the 1 :1 mass ratio to maximize the amount of NdFeB content while preserving the ease of fabrication.

B. Surface properties of the robot

The native surface of the composite elastomeric material is hydrophobic and mi- croscopically rough. Static advancing (1 16° ± 3°) and receding (78° ± 2°) water contact angles were measured by the sessile droplet method through an automated goniometer routine (Kruss DSA100 and accompanying software) that continuously increases and decreases the water droplet volume, respectively. Surface roughness (R _a = 0.63 ± 0.02 μιτι and R _z = 4.37 ± 0.28 μιτι) was measured by laser interferometry (Keyence VK-X200). The as-fabricated robots' surface was not subjected to any treatment before the experiments.

C. Biocompatibility of the robots

While the Ecoflex-10 polymer matrix is itself biocompatible, the presence of NdFeB microparticles within the matrix makes the current composition of our soft ro- bots only partially biocompatible. Full biocompatibility of our soft robots can be recovered by sealing the microparticles safely within the robots. For example, this can be achieved by coating and encasing the robots within an additional, outer thin shell of Ecoflex-10 polymer. This solution would also preserve the surface properties discussed above, as the native Ecoflex-10 is itself hydrophobic and has a microscopically rough surface. Furthermore, because unadulterated Ecoflex-10 has lower E than our composites shown in Fig. S1 d, the robot should still be sufficiently compliant to execute all of the proposed locomotion modes and functionalities as long as the additional layer is sufficiently thin. To ensure that we can predict the behavior of such robots, it would then be necessary to characterize the Young's modulus of this new composi- tion.

D. Residual strain energy

We hypothesize that there are two lowest energy configurations because our fabrication process is not perfect and the pre-stress induced from the demolding process will make the robot retain some strain energy. One consequent effect of this phenom- enon is that the proposed robots will exhibit two different rest state curvatures when no B is applied (Fig. S1 e). The first rest state curvature, which is relatively smaller, is illustrated in Fig. S1 e(l, V) while the second rest state curvature, which is larger, is shown in Fig. S1 e(lll). We found that the robot can alternate between these two rest state curvatures when it undergoes different shape-change. In particular, the relatively flat curvature can be restored after the robot recovers from the shape specified in Fig. S1 e(IV), and the larger residual curvature rest state is activated after the robot recovers from the shape specified in Fig.S1 e(ll). Our experiments show that these two different rest state curvatures can be switched repeatedly.

Having two rest state curvatures implies that the robot has two distinct lowest energy configurations. By inducing a shape-change to alter the amount of strain energy within the robot, the robot can switch between these two lowest energy configurations after B is removed. Theoretically, if there is no residual strain energy, the lowest energy configuration of the robot should be perfectly flat.

As mentioned in the foregoing and the captions in Fig. 19, the residual strain energy would also induce small errors in the predicted robot shapes. Since this fabrication uncertainty may vary across different robots, it is expected that different robots would have slightly different errors in their predicted shapes. To reduce such uncertainties, we always use the same batch of robots to characterize each locomotion mode (SI sections S4-10).

Despite having residual strain energy, we remark that our robot is still able to achieve multiple modes of locomotion. This is supported by our multi-locomotion demonstration in Fig. 3 where we show that one single robot is able to realize all the proposed modes of locomotion. In the future, we will account for these effects via cali- brating B and we will also improve our fabrication process such that the fabricated robot will have less residual strain energy.

S2 - Magnetic actuation setup

A. Experiment setup

The magnetic actuation setup is composed of three orthogonal pairs of custom- made electromagnets, and it has an inner chamber size of 65 mm χ 65 mm χ 65 mm (6). The input currents driving the electromagnets are specified by software control signals through a custom electronic board. We use a calibrated magnetic actuation matrix to map the control current input to the actual magnetic field B in the working space. With respect to the center of working space, which coincides with the center of the reference system defined in Fig. 6, the 95% homogeneous field region is 16 mm wide along the X- and Y-axis and 50 mm wide along the Z- axis. All demonstrations of individual locomotion modes were conducted in this homogenous field region, except for the multimodal locomotion experiments in Fig. 3a, d, e, for which the maximum dis- placement of the robot was around 30 mm. In these situations, we would recalibrate our actuation matrix to produce a relatively spatially uniform B for our robot. All characterizations in sections S4-S1 1 were conducted in the homogenous field region, except for the characterization of the jumping and jellyfish-like swimming, in which a pair of Helmholtz coil is used for simplicity.

Despite our efforts to provide spatially uniform fields, uncontrolled spatial- gradients in these regions may still exist. For example, we have observed unwanted magnetic forces pulling the walking robot towards the center of the workspace (Fig. 3d). While these forces have induced unwanted disturbances on the walking locomo- tion, we remark that the robot is able to walk successfully even against the gradient (see SI section S13 for more discussion). We further remark that the current setup was designed for proof-of-concept purposes, showing the possibility of creating highly mobile miniature soft robots that have multiple modes of locomotion. Electromagnetic coil setups with larger uniform workspaces can well be conceived to address these issues and render the whole robotic system more suitable for potential biomedical applications described in SI section S14.

B. Coordinate frames

Here the global frame is defined by the magnetic actuation setup's coordinate frame, and we shall define the axes in this frame with capital letters, e.g., X-, Y- and Z- axes, and XY plane. To distinguish this global frame from the robot's un-deformed local frame shown in Fig. 19, we assign the axes of the local frame with lower case letters, e.g., X-, y- and z-axes, and xy plane. As the soft robots can be steered by B, we can orientate the robot's local frame with respect to the global frame accordingly (see SI section S3B). For example, we can make the robot's local xy plane parallel with the global XY plane by applying a B that only has components in the XY plane.

Except for the main text and the shape change analyses in SI section S3A, which are more intuitive to be described in the robot's local frame, all other analyses and characterizations in the SI are described in the global frame (SI sections S4-1 1 ). In these analyses and characterizations, we have orientated the robot's local xy plane to be parallel with the global XY plane.

S3 - General theory

Here we provide the general theory, which describes how our soft robots respond to external magnetic fields. We begin by providing a quasi-static mathematical analysis of how the actuating fields can generate desired time-varying shapes for a magnetically actuated soft robot (section 3A). Subsequently, we describe how the actuating fields can be used to rotate the body of the robot (section S3B), and also how to use our theory to determine the operating range of B (section S3C). SI sections S4-1 1 then propose more detailed models and discussions for each locomotion mode. A scaling analyses summary for the soft robots is finally presented in SI section S12 and Table S4. We also provide a nomenclature for all the mathematical symbols in Table S1 and a summary of our robot's physical properties in Table S2.

We remark that the analyses and discussions presented here and in the following SI sections aim to facilitate future studies on soft-bodied locomotion. The analyses can also provide design guidelines, which can be used as the first step towards designing and optimizing miniature soft magnetic robots that have multi-locomotion capabilities.

A. Quasi-static analysis

Given the magnetization profile m prescribed in Fig. S1 b, here we use quasi-static analysis to describe how we can use the external magnetic fields,

to control our robot's shapes. The magnetization profile, m, along the robot's length, s, is described as:

The variables m _x and m _y represent the components of m along the x- and y-axis, respectively, while m and <¾ represent the magnitude and spatial angular frequency of m(s), respectively, with Using the above m(s) and without any loss of gener

ality, we use the robot's local reference frame in Fig. S1 b to describe B, which is a function of time, t. At a given t, the interaction between B and m(s) creates a spatially varying stress in the soft robot to create desirable deflections (Fig. S3). These static deflections on an element can be described by the Euler-Bernoulli equation for beams, i.e., the bending moment where the rotational deflection along the robot,

6(s), is expressed explicitly with respect to actuating magnetic torque, _m , as ³ _: where

The variables = hw, I = h ³ w /12 and E represent the robot's cross-sectional area, second moment of area and Young's modulus, respectively, while R is the standard z- axis rotational matrix that accounts for the change in direction of m due to the robot's rotational deflection:

Physically, the left side of Eq. (S3.2) represents the magnetic actuation that deforms the robot, while its right side is a function of the resultant deflections. Since we do not apply magnetic forces onto the robot, the shear forces along the robot are null at steady-state and are therefore neglected for the quasi-static analyses (Fig. S3). For the locomotion modes that have small deflections, such as the undulating swimming, Eq. (S3.2) can be simplified by approximating R with a 3 3 identity matrix, s (along the body frame's x-axis), and where y refers to the vertical deflection.

With these simplifications, Eq. (S3.2) then becomes:

Using Eq. (S3.5), the robot's shape under small deflection condition can be de- scribed via the linear combination of two principal shapes, resembling respectively a sine and a cosine function (Fig. 19b(ll,lll)). The first principal shape can be enforced when the direction of B is parallel to i.e., B = Then, Eq. (S3.5) can be triple integrated with respect to x to become ⁴ :

where B represents the magnitude of B, and / i, k ₂ and / ₃ are the constants of integration, valid for all B values as long as the robot's deflections are small in

cluding the case of B = 0. When B = 0, the LHS of Eq. (S3.6) has to be zero for all x, because the robot does not deform in the absence of B, and this implies that Eq. (S3.6) becomes:

As x can vary from 0 to L, Eq. (S3.7) can only be satisfied across all x if and only if / i = k ₂ = k ₃ = 0. Therefore, Eq. (S3.6) can be rewritten as:

Similarly, when B is parallel to i.e., B = Eq. (S3.5) can be triple integrated to become:

and using the same reasoning as for / i, / ₂ and / ₃ , we deduce that Eq. (S3.9) can only be valid across all x if and only if k ₄ = k ₅ = k ₆ = 0. Therefore, Eq. (S3.9) becomes:

As mentioned in the main text, if the direction of B is not aligned along the principal directions, the robot will assume a shape that can be described by the weighted superposition of Eqs. (S3.8) and (S3.10). While Eqs. (S3.8) and (S3.10) can be used to fully describe the small deflection shapes used in undulating swimming and menis- cus climbing, they are insufficient to describe the large deflection shape change necessary for rolling, walking, jumping and jellyfish-like swimming. The large deflection shapes are achieved by applying a large magnitude B parallel or anti-parallel to To predict the robot's shapes under such conditions, we have to solve Eq. (S3.2) numerically with the following free-free boundary conditions:

Physically, these boundary conditions imply that there is zero bending moment applied on the free-ends of the beam. The numerical solutions reveal that the robot creates a large deflection configuration that resembles a 'C'-shape when we apply a large magnitude B parallel (Fig. 19b(IV)). On the other hand, when we use the same boundary conditions and apply a large magnitude B parallel to the numerical solution of Eq. (S3.2) yields a 'V'-shape (Fig. 19b(V)). In both configurations, the curvature increases as B is increased. Further- more, if the direction of B is not aligned along the principal axis shown in Fig. 19b(IV,V), the robot will have a rigid-body rotation until its M _ne t is aligned with B. At the end of this rotation, the robot will assume its 'C- or 'V'-shape (see the next sub-section for more details).

B. Rigid-body rotation

As the net magnetic moment of a robot tends to align with the external magnetic fields, it is possible to use the direction of the fields to control the orientation of the magnetic robot ^{1 ,5} . While our robot does not have a net magnetic moment in the un- deformed state, a net magnetic moment arises once the robot deforms. Mathematically, the net magnetic moment, M _ne t, is expressed as:

Once equipped with M _ne t, we can create rigid-body torques on the robot by controlling B. This magnetic torque, XRB, can be represented as:

Physically, Eq. (S3.13) implies that as long as M _ne t is not aligned with B, a magnetic torque will be applied onto the robot. If there are no other external torques acting on the robot, the magnetic torque will eventually align M _ne t along the direction of the applied B at steady state conditions. In the following sub-sections, we will discuss how we can rotate the robot about its z- and y-axes.

I. Rotation about z-axis

Here, we will explain why the robot is unable to rotate about its z-axis when it has small deflections (i.e., B is small) but will be able to do so under large deflection condi- tions (i.e., B is large). We begin this discussion by showing that the robot's M _ne t will always be aligned with B under small deflection conditions, and this in turn will conclude that there will be no rigid-body magnetic torques acting on the robot, making the robot unable to rotate about its z-axis. To establish this mathematical proof, we first show that when the robot is subjected to a B that is parallel to the principal directions shown in Fig. 19b(ll-lll), its resultant deformation will generate an M _ne t parallel to the corresponding principal directions. Subsequently, by using the superposition principle, we can generalize and conclude that the robot's M _ne t will always be aligned with the applied B under small deflection conditions.

Based on Eq. (S3.8), the robot will produce a cosine shape when B is parallel to its first principle direction, i.e., Since the deformation of this shape is small, its rotational deflection can be approximated as:

Using the small rotational deflection representation in Eq. (S3.14), we can approximate cos0 and sin# to be 1 and Θ, respectively. These approximations will in turn simplify the resultant rotational matrix R in Eq. (S3.12) and Eq. (S3.4) to:

By substituting Eq. (S3.15) into Eq. (S3.12), the deformed robot's M _ne t can be expressed as:

Based on simple trigonometry identities, can be re

written as:

By substituting Eq. (S3.17A-B) into Eq. (S3.16) and rearranging the equation, the deformed robot's Mnet can be expressed as:

As the first component on the right side of Eq. (S3.18) integrates the harmonic functions over one full wavelength, this component is computed to be zero. Therefore, by completing the integration for the second component on the right side of Eq. (S3.18), this equation can be simplified into:

We remark that the direction of M _ne t in Eq. (S3.19) is parallel to the first principal direction of

Next, we repeat similar derivations when B is parallel to the second principal direction, i.e., Based on Eq. (S3.10), the robot will produce a sine shape with this B. Assuming small deflection conditions, its rotational deflection can be approximated as:

and the corresponding rotational matrix in Eq. (S3.12) becomes:

By substituting Eq. (S3.21 ) into Eq. (S3.12) and use steps similar to those presented in Eq. (S3.16-S3.18), the deformed robot's M _ne t can then be expressed as:

and this shows that the generated M _ne t is parallel to the second principal direction of B. Hence, Eqs. (S3.19) and (S3.22) conclude that the deformed robot will not be able to rotate about its z-axis when a B with small magnitude is applied parallel to either of the principal directions, i.e., there is no rotation when the robot deflections are small. Because superposition principle is applicable for the analysis shown in Eq. (S3.14- S3.22), we conclude that under small deflection conditions (i.e., B is small), the robot will not be able to rotate about its z-axis. Hence in such situations, the robot will only experience a shape-change in the xy plane as described in the main text and SI section S3A.

Conversely, because the superposition principle will no longer be applicable for a large B, the robot's M _ne t in such situations will not necessarily be parallel with the ap- plied B. This eventually will make the robot rotate until its M _ne t is aligned with B, and the robot will assume either its 'C- or 'V'-shape configuration when it stops rotating. The 'C- or 'V'-shapes are the steady-state configurations because their generated Mnet is always naturally aligned with the applied B (as shown in Fig. 19b (IV-V)). Hence for large deflection cases, we can adjust the direction of B in the XY plane to rotate the 'C- or 'V'-shape robot about its z-axis (Fig. 19c). The ability to control the robot's angular displacement about its z-axis is essential for achieving several modes of locomotion, as it allows the robot to produce the necessary tilting for jumping and walking, as well as enabling the rolling locomotion (SI section S5).

Finally, we will also discuss what happens when a mid-range B is applied, i.e., a B that is neither sufficiently small nor big. In this situation, two possible scenarios could occur. When a mid-range B is lower than a specific threshold, the deformed robot will not be able to generate an effective rigid-body rotation. As a result, the obtained shape of the robot can be predicted by solving Eq. (S3.2) numerically. On the other hand, if the mid-range B exceeds this threshold, the robot would rotate until its M _ne t is aligned with B and it will assume either its 'C or 'V'-shape configuration in steady-state conditions. Theoretically, this threshold represents the minimum required B for the deformed robot to generate a M _ne t that is non-parallel to the applied field. However, in practice as our fabrication process is not perfect, it is not easy to quantify this threshold accurately. Hence, avoiding mid-range B and employing either small or large B would be the best to ensure that the robot is able to achieve all modes of locomotion. Indeed, we show that we are able to perform all modes of locomotion by using ranges of B that are sufficiently small or large.

II. Steering strategy (rotation about y-axis)

As the robot's M _ne t always resides in its xy plane (Fig. 19b), this implies that when a B that only has XY plane components is applied on the robot, a magnetic torque can re-orientate the robot until its local xy plane becomes parallel to the global XY plane shown in Fig. 6. Similarly, after making the robot's xy plane to be parallel with the global XY plane, we can use B _z to create a magnetic torque that can rotate the robot about its y-axis and steer it along other desired directions. In contrast to rotating the robot about its z-axis, we can always rotate the robot about its y-axis regardless of the magnitude of B. More information regarding how to use external magnetic fields to steer a magnetic robot with a net magnetic moment along a desired direction can be found in references 1 and 3. C. Predicting the operating range of B

The theoretical models presented in section 3A-B can help to determine the operating ranges of B such that the robot can generate the required time-varying shapes to enable different modes of locomotion. For locomotion modes enabled by small B, we can use Eqs. (S3.8) and (S3.10) to predict the required B for generating the necessary time-varying shapes. However, a recursive approach has to be used for determining the required B for locomotion modes that require a large B. This recursive approach is implemented by first specifying a B, and substitute it into Eq. (S3.2) to numerically solve for the curvature of the 'C- or 'V'-shape configuration. If the obtained curvature is unsatisfactory, we will change another B, and repeat the same procedures until a suitable B is obtained.

S4 - Analysis for jumping

To analyze the jumping locomotion, we first specify the required B signals (section S4A). Subsequently, we provide a theoretical model, which can predict how high the robot can jump when we vary B (section S4B). Using this model, we then discuss how the jumping height of the soft robot is affected by its dimensions (section S4C).

A. Experimental details

For the robot to jump along a desired direction (Fig. 20h, also in Fig. S4a), we specify the following B sequence (Fig. S5a):

1 . The robot is first curled into a 'C'-shape configuration by specifying B = 16.5 mT and a = -90° (0 ms in Fig. S4a).

2. Subsequently, we use a step function to change the direction of B, rotating its direction 162° clockwise. This step function simultaneously induces a rigid-body rotation (clockwise in this case) as the robot's M _ne t tends to align with the direction of B, as well as the robot's shape changes shown in 0 - 13.5 ms in Fig. S4a.

3. We maintain B (B = 16.5 mT, a = 72°) until the ends of the robot hit the ground (13.5 ms in Fig. S4a), allowing the robot to jump (13.5 - 16.9 ms in Fig. S4a).

As mentioned, such B sequence prompts two different mechanisms as the robot performs directional jumping. The first mechanism is the rigid-body rotation of the ro- bot, induced by the rigid-body torque created from the interaction of M _ne t and the flipped B (0 - 13.5 ms in Fig. S4a). The second mechanism is the shape-change, as the robot changes from its initial 'C'-shape (0 ms in Fig. S4a) to a flatter sheet (10.2 ms Fig. S4a) and then back to its 'C'-shaped configuration (10.2 - 13.5 ms in Fig. S4a). Both mechanisms help to increase the momentum of the robot until the robot strikes the substrate. Upon impact, the substrate then induces both a vertical and horizontal momentum for the robot to jump. In the following discussion, we focus on the shape-change mechanism, for which we have devised a 'straight' (vertical) jumping strategy (Fig. S4b). Using robots with /3R = - 90°, the required setup for experimental characterization reduces to the single pair of electromagnetic coils along the Y axis (Fig. 6). Robots with /3R = 45° can also per- form straight jumping, though they require a more complex experimental setup, i.e., two pairs of electromagnetic coils to generate a step B along the a = 135° direction, and this may induce unnecessary experimental errors. To further ensure a spatially uniform B, we adjusted the distance of the two coils to form a Helmholtz configuration. Using these experimental conditions, the robot jumps when we specify a step change of B along the Y-axis of the global reference system. Mathematically, B can then be simply expressed as

B. Theoretical model

Here we provide a model for the straight jumping strategy. This model is based on an energy conservation approach, which first analyzes the energy states of the soft robot at two salient instants: 1 ) when it is initially flat on the substrate, and 2) the moment when it leaves the substrate. Based on the principle of energy conservation, the soft robot gains energy from the work done by the magnetic torque, W, between these two states. This energy is then re-distributed into three components: the change in strain energy, AS, and kinetic energy, AK, for the robot, and the frictional losses, f , during the jump of the robot. Mathematically, this implies that

The frictional losses include energy loss from contact friction with the substrate, air drag, and viscoelastic losses within the soft robot. To simplify our analysis, here we will neglect the effects of frictional losses and assume that f _L is equal to zero. Equation (S4.1 ) can therefore be simplified into

Based on the Euler-Bernoulli equation, AS can be expressed as

where 0 (s) represents the final angular deformation of the robot before it leaves the substrate and we assume that the initial shape of the robot is flat, i.e., 0 _£ (s) = 0.

On the other hand, the magnetic work done 1/ can be expressed as a double integral:

The inner integral represents the work done by the magnetic torque on an arbitrary infinitesimal volumetric element along the robot (Fig. S6). The upper limit of this integral represents the final angular deformation of the element before the robot leaves the substrate and the lower limit 0 _£ (s) = 0. Therefore, 1/ can be obtained by the outer integral as we sum up the total work done over all the infinitesimal elements along the robot. By using Eqs. (S3.3-S3.4), the inner integral can be evaluated to be

Because the magnetization profile components, m _x and m _y , and the actuating fields, B _x and By are independent of Θ, Eq. (S4.5) can be evaluated to be

Finally, by substituting Eq. (S4.6) and our magnetic actuation signals

= [0 B 0] ^T ) back into Eq, (S4.4), we can express W as

Physically, Eq. (S4.7) implies that the work done is proportional to the magnitude of the magnetization profile and the actuating field, and it is also dependent on the final shape of the robot before it leaves the substrate. By substituting Eqs. (S4.7) and (S4.3) back into Eq. (S4.2), we obtain a mathematical expression for ΔΚ. Once ΔΚ is known, we apply the principle of energy conservation again, this time from the moment when the robot leaves the substrate till the moment when it reaches its maximum height, H _max . At the maximum height, we assume all the kinetic energy used for jumping is fully converted into gravitational potential energy such that

where M _Ro b and g represent the mass of the robot and the gravitational constant, respectively. We note that Eq. (S4.8) is a simplified model because it does not account for frictional losses and also not all the kinetic energy in Eq. (4.2) is used for jumping, i.e., some of the kinetic energy is converted to excite other modes of vibration of the robot. By substituting Eqs. (S4.3) and (S4.7) into Eq. (S4.8) and by rearranging the equation, H _max can be obtained as

Using Eq. (S4.9), we can then predict the theoretical H _max of the soft robot across different B. Equation (S4.9) shows that H _max \s independent of w but its highly non-linear nature makes it difficult to derive a simple scaling law for h and L. Despite so, we will show in the following sub-section that this model is able to predict that the sample robots in Fig. S7 can jump higher if they have a smaller h (Fig. S7a, b) and a larger L (Fig. S7a.c).

C. Discussion

To investigate how the jumping performance varies when the dimensions of the soft robot are changed, we fabricated six additional robots, each having a different L, w or h from the original robot shown in Fig. 19a. For these additional robots, we have measured their maximum jumping height for 13.5 mT < B < 23.5 mT. We used image detection techniques in MATLAB and standard spline fitting methods to extract from the experiments the final shape of the robots before they left the substrate, i.e., the Of of the robots. Experimental and model-predicted values of H _max are presented in Fig. S7. Because Of may vary in each experimental point, the model predicted points in Fig. S7 also have error bars.

Generally, while the model in Eq. (S4.9) agrees with the experimental trends, the predicted values of H _max in Fig. S7 overestimate the actual values as they assume that there are no frictional losses and all kinetic energy of the robot is used for jumping.

The experimental dependency of jumping height on B parameterized by the dimensions of the robots is shown in Fig. S8. Based on the experimental data, robots can jump higher if they have a longer L or a thinner h and there are no observable effects when w is changed. This trend agrees with the model specified by Eq. (S4.9) (see model predicted curve in Fig. S7).

In summary, our experimental results suggest that the soft robot can jump higher if we increase L or decrease h for the straight jumping strategy. The predictions of the theoretical model agree with the experimental data.

S5 - Analysis for rolling

Here we analyze the rolling locomotion by first specifying the required B sequence (section S5A) and then providing a theoretical model, which can predict how fast the robot can roll as we vary the rolling frequency f of B (section S5B). Finally, we discuss how the rolling speed Vroii of the robot is affected when we vary the robot's dimensions (section S5C).

A. Experimental details

A schematic drawing for the robot rolling in Fig. 20e is illustrated in Fig. S10. This locomotion is realized by using the following B sequence in the XY plane (Fig. S9a):

1 . An initial B (B = 18.5 mT, a = 0°) is imposed to curl the robot into a 'C'-shape as shown in Fig. S9b.

2. The deformed robot possesses M _ne t as specified by Eq. (S3.12), allowing it to be steered by a rotating B in the XY plane. The rotating frequency is 5 Hz in Fig. 20e.

3. As B starts to rotate clockwise, a rigid-body magnetic torque, TR ₀ H, is induced onto the robot as the robot's M _ne t tends to align with B. This torque makes the robot roll.

4. Once B stops rotating, the robot's M _ne t fully aligns with B and TR ₀ H is reduced to a null vector. Hence, the robot will stop rolling.

Mathematically, the sequence of B described above can be represented as:

where f represents the frequency of the rotating B and it is ranged between 0 < / < 40 Hz in all of our experiments. To fairly evaluate V _r0 ii for robots with differing dimensions, we had adjusted B for each robot such that all of them could produce the same curvature. This adjustment is described in the next subsection.

B. Theoretical model

To enable the rolling locomotion, the magnitude of the rotating B has to be large to make the robot curve into a 'C'-shape. For the original robot shown in Fig. 19a, a B larger than 10 mT is selected to guarantee a sufficiently large robot deflection.

We assume the rolling locomotion akin to that of a wheel rolling on a flat sub- strate. This assumption is supported by the observation that the robot can curl into a 'C'-shape that resembles a circular wheel when we apply the rotating B. We have used MATLAB 2016 (boundary value problem solver) to numerically solve the governing Eq. (S3.2), and we then computed the net magnetic moment, M _ne t, of the robot in the 'C'-shape configuration with Eq. (S3.12). In the rolling locomotion, M _ne t becomes a function of the rotational displacement of the robot, _R (Fig. S10), and mathematically it can be expressed as

where the scalar M _net represents the magnitude of the net magnetic moment.

Once B starts rotating, a rigid-body torque, τ& ₀ ιι, is applied on the robot as its net magnetic moment tries to align with the field. Mathematically, TR ₀ H can be written as

The magnitude of TR ₀ // depends on the M _net and B, and on the angle between them. As the robot rolls, the ground induces a pushing force, F, which propels the translation of the robot (Fig. S10). The translational equation of motion of the robot is

where Cr represents the translational damping. Likewise, the equation of motion of the robot about the z-axis can be expressed as

where CR, r _eff , and J represent the rotational damping coefficient, radius of wheel and moment of inertia of wheel, respectively. Assuming that the robot does not slip as it rolls on the substrate, the kinematic relationship between the x and φ _κ can be expressed as

By substituting this no-slip condition into Eq. (S5.4), multiplying it by r _eff , and finally adding it to Eq. (S5.5), the governing dynamic equation for the rolling locomotion can be expressed as

where Since the size of the robot is relatively small, its small / _eff enables the robot to quickly reach a steady-state velocity, which can be rewritten as

A common phenomenon for miniature, magnetically actuated rotating robots is that the angle between the actuating field and the robots' angular displacement be- comes constant at steady-state speed ^6,7 . Mathematically, this means that sin (2nft - φ _κ ) becomes a constant value in Eq. (S5.8). Equation (S5.8) also reveals that the maximum angular velocity of the robot occurs when the angle between the field and the robot's angular displacement becomes 90°, i.e., since sin(27i: t—

0 _R ) = 1. This maximum angular velocity is known as the step-out frequency for the soft robot, and it is represented by

Once the rotating frequency of B exceeds the step-out frequency, i.e.,

the frequency of B is no longer synchronized with the robot's angular veloc¬

ity and it is difficult to predict φ _κ . The discussion for such conditions will, however, be beyond the scope of this paper.

We experimentally validated that our robots do not exceed the step-out frequency under the control signals we used. Therefore 2nft - φ _κ = k, where k is the constant lag angle between B and φ _κ , and it does not change with time once the robot reaches its steady-state speed. By differentiating this expression in terms of time, this implies that

Based on Eq. (S5.10) and the no-slip condition in Eq. (S5.6), the governing equation for our rolling locomotion becomes

r _eff can be determined by first solving Eq. (S3.2) numerically using MATLAB 2016 (boundary value problem solver) to determine the rolling shape, i.e., the 'C'-shape, and then determine the parameters of the circle that best fits it. We remark that the effec- tive circumference of the fitting circle, i.e., is generally approximately 10% higher than the length of the robot. A simple scaling analysis can also be derived based on

Eq. (S5.1 1 ) if we approximate and this will suggest that

As m, B, E and other constants are independent of the dimensions of the robot, we can conclude that a higher B is required to roll stiffer robots, i.e., robots with smaller L or larger h ³ . Indeed, to enable a fair comparison of rolling locomotion performance for robots of different dimensions, we have adjusted B in such a way that all of the robots could produce the same curvature. For the robots with length L _new different from the reference L, the same curvature could be obtained when the magnitude of B, became ³ :

Likewise, we could obtain the same curvature for robots that had thickness h _new by changing the magnitude of B, to

Finally, varying w does not change the curvature if B is unchanged. Therefore, we used Eq. (S5.12-S5.13) in our experiments to adjust the corresponding B for robots with differing dimensions, allowing all of the robots to create a curvature equivalent to that of the original robot under B = 18.5 mT.

C. Discussion

Here we discuss how the rolling speed V _r0 ii varies with the frequency f of the rotating B and the dimensions of the robot. For this purpose, we fabricated six additional robots where each of them has either a different L, w or h compared to the original robot.

The experimental data presented in Fig. S1 1 show that the V _r0 ii has a linear rela- tionship with f at low values, e.g., f < 10 Hz in Fig. S1 1 a-b. Generally, our theoretical model agrees with the data, though it underestimates V _r0 ii. This discrepancy could be justified by the fact that the actual 'C'-shape configuration of the robot is not a perfect circle, as assumed in the model. Moreover, the standard deviation of V _r0 ii increases when f becomes high, e.g., when 20 Hz < / < 40 Hz in Fig. S1 1 a-b. This occurs be- cause for such high f the ends of the robot tend to impact the substrate with such a large force that it makes the robot jump. To avoid this less predictable behavior, we recommend constraining within its linear regime.

The experimental data for the dependency of V _r0 ii on f parameterized by the dimensions of the robots, presented in Fig. S12, show no significant variations in versus w and h. The V _r0 ii, however, increases when the robots have larger L. This trend generally agrees with the predictions of Eq. (S5.1 1 ), as r _e ff becomes bigger when L increases.

In summary, the rolling speeds of the robots have a linear relationship with a low f. Robots with a larger L can roll faster. Rolling is also much more controllable for rela- tively low f {e.g., f ≤ 4 Hz for our robots). Finally, according to Eqs. (S5.12) and (S5.13), it is easier to induce the curvature necessary for rolling (i.e., we can induce the same curvature using a lower B) for robots with larger L or smaller h. In this respect, we can use the same B to achieve the same curvature for robots that have different w. S6 - Analysis for walking

In this section, we first specify the sequence of B required for walking (section S6A) and then provide a model that can predict the walking speed of our robot as we vary f and B _max , i.e., the maximum B value in one walking cycle (section S6B). Finally, we discuss how the robot dimensions affect its walking speed in section S6C.

A. Experimental details

Driven by the B sequence shown in Fig. S13a in the XY plane, the soft robot can walk rightwards as shown in Fig. 20f. The B sequence can be described as follows:

1 . The robot first anchors on its front end and starts to tilt forward. As it rotates, it also starts to deform into a 'C'-shape (0 ms - 62 ms in Fig. 20f). To realize this motion sequence, the direction of B is fixed at a = 150° while B is increased linearly until it peaks at B _max = 10 mT.

2. The robot maintains the 'C'-shape and begins to anchor on its back end such that it can tilt backwards (62 ms - 122 ms in Fig. 20f). To realize this motion, the direction of B changes linearly from a = 150° to a = 120° while maintaining B at 10 mT.

3. Finally, the robot extends its front end, allowing it to make a net stride forward (122 ms - 210 ms in Fig. 20f). To realize this motion, the direction B is fixed at a = 120° while its magnitude is decreased linearly from 10 mT to 0 mT.

To investigate how the walking locomotion is affected by the dimensions of the soft robot, we fabricated six additional robots that have either a different L, w or h compared to the original robot. Using these robots, we conducted two types of experiments, where we respectively evaluated the walking speed of the robots, V _wa ik, versus different f (from 2 Hz to 20 Hz) with B _max fixed at 10 mT, and V _wa ik versus different B _max (from 8 mT to 16 mT) with f fixed at 5 Hz.

B. Theoretical model

Based on our walking strategy, the net stride of the robot in one walking cycle,

Swaik, is a function of the difference between the robot's curvature in the states shown at times 0 ms and 102 ms in Fig. 20f. More explicitly, S _wa ik can be expressed as

where Si and S ₂ represent the distance between the ends of the robot as illustrated in Fig. S14. Because B = 0 in Fig. S14a, Si can be approximated as L. On the other hand, S ₂ in Fig. S14b can be computed by first solving Eq. (S3.2) numerically with MATLAB 2016 (boundary value problem solver). We assume that the soft robot can respond to the actuating fields when the frequency of B is low. In these regime, the walking speed, V _wa ik, can be approximated as the product of S _wa ik and the walking frequency f

Validity of Eq. (S6.2) assumes no-slip conditions at the anchored ends. Equations (S6.2) and (S3.2) suggest that V _wa ^ _k does not scale with w, and they also suggest that ^v wa k ^{nas a} non-linear relationship with L and h, as S2 can only be solved numerically by Eq. (S3.2). However, Eq. (S3.2) does suggest that robots with larger L and thinner h can deform more to create a smaller S2. This in turn suggests that Eq. (S6.2) predicts that robots with larger L and thinner h will walk faster.

C. Discussion

Figure S15 compares the experimental data with corresponding theoretical predictions, based on Eq. (S6.2) and Eq. (S3.2), for V _wa i _k versus f. The experimental results suggest that V _wa i _k has a linear relationship with f at low frequencies. However, for higher values of f, there appear instances where V _wa ik has a non-monotonic trend, e.g., Vwaik drops to zero after 14 Hz in the L = 333 μιτι case (Fig. S15b(ll)). This happens because, under such conditions, the ends of the robot approach the substrate so quickly (at the end of Step 2 in SI section S6A) that upon contact the robot is bounced backwards. Such bouncing reduces the effective net stride per walking cycle, making the robots slower, and in the case shown in Fig. S15b(ll), waik becomes 0 mm/s.

Figure S16 compares the experimental data with corresponding theoretical predictions, based on Eq. (S6.2) and Eq. (S3.2), for The experimental results show that Vwaik correlates positively with S _max . This agrees with our theoretical model, which predicts that a larger S _max can induce a larger net stride, i.e., (L -S2), in each walking cycle.

The experimental data for the dependency of V _wa ik on f and S _max parameterized by the dimensions of the robots are respectively shown in Fig. S17 and Fig. S18. The results in Fig. S17 suggest that, in the linear regime where f < 1 1 Hz, more compliant robots (i.e., with larger L or smaller h) walk faster when S _max is fixed. This happens because the more compliant robots can produce a larger net stride per walking cycle. Additionally, the robots with larger L are expected to walk faster thanks to a relatively larger body-length. Our experimental data do not show any obvious correlation between Vwaik with w.

In general, while the theoretical predictions in Eq. (S6.2) and Eq. (S3.2) tend to overestimate V _wa ik, they agree with the trend of the experimental data (Fig. S15 and S16). Indeed, our model predicts that robots with smaller h and larger L walk faster while w does not affect V _wa ik- The predicted speeds are larger than the experimental speeds because the robots are observed to slip at their anchor points on the substrate, and slippage reduces the net stride of the robots in each walking cycle. Furthermore, due to fabrication uncertainties (section S1 D), the robots are not perfectly flat even for B = 0. Therefore, also by approximating Si with L our model tends to overestimate the net stride in each walking cycle.

S7 - Analysis for meniscus climbing

In this section, we analyze the robot's meniscus climbing by specifying the re- quired sequences of B for this locomotion (section S7A) and then proposing a quasi- static, force-based model which can predict the minimum required field magnitude, Smin, to effectively let the robot climb up to the top of a positive liquid meniscus (section S7B). Finally, we discuss how the robot dimensions affect meniscus climbing in section S7C.

A. Experimental details

Driven by the B sequence of Fig. S19 in the XY plane, the water meniscus climbing shown in Fig. 20b that prompts a subsequent landing onto the adjacent platform, is achieved by the following steps:

1 . The initial B (B =18.5 mT, a = 315°) is imposed to let the robot assume a 'C'-shape with upward curvature. The robot achieves in 520 ms a stable position along the meniscus.

2. B is subsequently rotated counter-clockwise toward a = 157° to make the robot climb further up. After 4.6 s of B is on, the robot reaches contact with the wall of the platform. The subsequent rotation progressively pushes the upper portion of the robot out of water while being supported by the curled lower portion of the robot on the meniscus (4.8 s - 12.1 s).

3. Finally, B is turned off and the robot rests on the edge of the adjacent platform.

The above description uses a 'C'-shape configuration to climb the meniscus. If there is no need to push the upper portion of the robot out of the water, we remark that the meniscus climbing locomotion can also be realized when the robot has cosine- shaped small deflection (Fig. 19blll). In the following section, we propose a quasi- static model of the soft robot's water meniscus climbing, which describes the process till the robot touches the wall of the platform. The corresponding experiment we performed is shown in Fig. S20. A B pointing upwards and of increasing magnitude was imposed to progressively curve the robot, and the smallest B magnitude S _m in which made the robot touch the wall at the top of the meniscus was recorded. As the robot's M _ne t is always aligned with B under small deflection conditions (section S3B), the following model does not consider magnetic torque.

B. Theoretical model

When contacting water from air, our soft robots get almost completely immersed in water (Fig. S20) and strongly pinned at the water/air interface, as confirmed by the conformal attachment of the triple contact line to the robot's perimeter during meniscus climbing (see e.g., Fig. 3d).

The robots stably reside at the water/air interface thanks to the generalized Archimedes principle, i.e., the combined effect of buoyancy, acting on the surface of the robots in contact with water, and of surface curvature, acting on their contact perimeter. The generalized buoyancy force is quantified by the total weight of water displaced by the robots, which includes the water displaced both by the robots' body (i.e., inside the contact line) and by the menisci along their perimeter (i.e., outside) ⁸ . Given the conformal contact line pinning along its perimeter, a floating body can move upward along a positive liquid meniscus provided that the quasi-static work the body makes on the system, by extending and/or curving the water interface, compensates the energy required to climb the gravitational potential.

Here we first describe a two-dimensional force-based analytical model for an arc- shaped, soft floating robot (arc-shaped model, in short) to explain the mechanics of meniscus climbing and predict the body curvature that allows the robot, with known geometrical and material parameters, to climb the height of a sloped water meniscus. We then derive a model for cosine-shaped robots (i.e., cosine-shaped model) to predict the minimal value of B, S _min , that makes the robot climb the complete height of the meniscus.

I. Arc-shaped model

The two-dimensional geometry of reference is sketched in Fig. S21 a. According to experiments, we assume the hydrophobic robot body fully immersed in water with only its top surface in contact with air. The triple contact line is consequently pinned at the perimeter of the top robot surface. Contact line pinning leaves the edge angles θ ₁ and θ ₂ undetermined, thus free to conform to boundary conditions within the limits imposed by canthotaxis, specifically by the Gibbs' criterion. This has two important consequences. First, edge angles are not imposed by material properties, and so they are largely independent of robot's pose (described by the tilt angle a) and curvature (by the radius of curvature R and the angle φ). As a corollary, the physical cause of menis- cus climbing is in our case not primarily related to surface tension effects. Second, by extending the canthotaxis sector, the robot's surface hydrophobicity enhances contact line pinning and thus the stability of the robot position on the water surface.

The mechanism underlying a soft robot's meniscus climbing is based on buoyancy. When assuming a curved shape, the robot (of density p _r ) displaces an additional volume of water (of density p _w ) from the surface (see Fig. 20b, reproduced in Fig. S19b). In presence of a positive water meniscus, the additional buoyancy enables the robot to move up along the interfacial profile. The five variables in the problem (e ₁ , e ₂ , a, R, and (p, see Fig. S21 a) can be determined by satisfying conditions for mechanical and hydrostatic equilibrium of the curved robot at a distance d from the adja- cent, vertical wall. The wall intersects the water surface profile, H(x), at a fixed contact angle θ ₀ . In experiments, both the vertical position H(0) of the contact line at the wall and θ ₀ were not affected by the presence of the robot (Fig. S20).

The shape of liquid interfaces in hydrostatic equilibrium is described by the Young- Laplace equation

The unperturbed profile H _u of the meniscus (i.e., in absence of the robot) is determined by the boundary conditions leading to

where we define the capillary length

The profile of the water surface H _r joining the wall to the closest robot edge is in turn determined by the Young-Laplace equation and the boundary conditions leading to

Assuming according to experiments, the following condition links d

to θ ₁ :

Mechanical equilibrium is imposed by normal and tangential force balances, with respect to the direction identified by α, and by torque balance. The three balances are respectively described by

where

and the distance of the centroid of the circular

segment from the center of curvature is In Eqs. (S7.5-S7.7), the first

terms represent the corresponding projection of the surface tension force, the latter term the effect of gravity and buoyancy. Finally, the conservation of robot length provides the last equation as

The system of nonlinear equations was solved numerically (MATLAB 2016, nonlinear system solver). In Fig. S22, we compare experimental and predicted values of robot's body curvature when first in contact with the wall at the top of the meniscus (i.e., for d→ 0). The actual shape of the deformed robot body was fitted with a circular pro- file in the same way as in SI section S5b.

II. Cosine-shaped model

To link body curvature to prescribed B, we derived a cosine-shaped version of the previous model to specifically describe a cosine shape robot lying at the top of the water meniscus, i.e., for d→ 0. The reference two-dimensional geometry is sketched in Fig. S21 b. In this case, we assume the ends of the robot pinned at the water-air interface, and that the latter has a constant (linear) slope that, from Eq. (S7.2), can be described as H'(0) = -cot # ₀ . From Eq. (7.4) we also get Under these conditions, horizontal force balance determines θ ₂ = θ ₀ . Accordingly, the tensile components of the surface tension acting on the ends of the robot are aligned with the water- air interface. Since the vertical projections of the surface tension force are of equal magnitude and opposite in direction, vertical force balance then imposes the equilibrium of weight and buoyancy of the robot. This implies the coincidence of center of gravity and center of buoyancy, as well as balance of torques. The buoyancy of the robot with cosine shape of magnitude D and spatial frequency <¾ can be calculated as

and the balance of robot buoyancy and robot's body weight leads to where the second factor is always positive for _r > p _w . Finally, by using Eq. (S3.8) the prescribed B _min to climb the entire meniscus height is computed as

The above equation shows that Comparison of experimental and model-predicted values of B _mm are presented in Fig. S23. C. Discussion

To discuss how the minimal prescribed B, B _mtn to climb the entire meniscus height, varies against the robot's dimensions, we include the experimental data of four additional robots that have either a different L, h compared to the original robot (Fig. 19a). As shown in Fig. S22, the arc-shaped model correctly predicts how L and h can impact the meniscus climbing performance: shorter and thicker robots require correspondingly smaller /? (i.e., larger curvature) to climb the full meniscus. Similarly, Fig. S23 shows that Eq. (S7.1 1 ) of the cosine-shaped model, specific for the robot in contact with the wall at the top of the water meniscus, can provide estimates of upper bounds on the values of B in accord with the experimentally recorded values.

With regard to the latter model, we remark that the simplifying assumption of constant linear slope of the water-air interface holds in principle for bodies of length L « L _c , for which the Taylor series expansion of the exponential meniscus profile (Eq. (S7.2)) can be approximated with its linear term. This condition is admittedly not quite met by our robots. In spite of this discrepancy, the model yields reasonable values for the prescribed B _m i _n .

In all cases, w is excluded from consideration, since the 2D nature of the models implicitly assumes infinitely wide robots. In experiments, the finite width of the actual robots introduces an additional curvature of the water surface, in the direction perpendicular to the 2D model plane (Fig. 3d). This transverse curvature is not contemplated in the model. Such inherent discrepancy may partly account for the discrepancy in the results. An additional approximation was introduced in the arc-shaped model by assuming that the additional volume of displaced water is equal to that of the circular segment subtended by the curved robot. This is admittedly a simplifying overestimate of the water volume additionally displaced from the surface, which may explain the model's systematic overestimate of R.

S8 - Analysis for undulating swimming

In this section, we analyze the undulating swimming by first providing the required sequence of B (section S8A) and then a model that can approximate our robots' swimming speed (section S8B). Finally, we discuss how the robot's dimensions can affect its swimming speed in section S8C. A. Experimental details

Similar to the rolling locomotion (SI section S5A), we use a rotating B in the XY plane to enable undulating swimming (Fig. S24). The key difference between these two modes of locomotion is that B has to be much smaller here {e.g., 1 mT to 5 mT) to prevent the robot from curling up into a 'C'-shape. With this rotating B, the robot is able to create an undulating traveling wave along its body such that it can produce an effective swimming gait similar to a Taylor swimming sheet. For the undulating swimming locomotion, the frequency of the rotating B, f, ranges from 40 Hz < f < 160 Hz, and B is expressed as

B. Theoretical model

The speed of the undulating swimming depends heavily on the highest rotational frequency of B that the robot can respond to. Therefore, we need to estimate the robot's mechanical bandwidth, i.e., its two lowest non-zero fundamental natural frequency ω _{η 2} and ω _{η 3} . This computation can be achieved by a simple vibrational analysis, which assumes small deformations within the robot. Using this assumption, we analyze an arbitrary infinitesimal element of the robot, dx, at time t (see Fig. S25). Based on Newton's force law, the equation of motion of such element along the vertical axis is given as

where C and v represent the damping coefficient and internal shear force, respective- ly, while p _r represents the robot's density and y refers to the vertical deflection. Physically, -Cwy and p _r Ay in Eq. (S8.2) represent the damping and inertial forces acting on the element, respectively.

In a similar manner, we can use r = Ja to obtain the rotational equation of motion for this element about its bending axis:

The variable M _b represents the bending moment acting on the element, and based on the Euler-Bernoulli equation, it can be approximated to be for small deflections.

Therefore, by substituting Eq. (S8.3) into Eq. (S8.2), we can obtain the governing equation for the beam:

We utilize the well-established mode-shape analysis to evaluate this partial differential equation. We first perform a free-vibration analysis to obtain the robot's mode shapes, i.e., we temporarily remove the damping forces and magnetic torques from Eq. (S8.4). Subsequently, we use the variable separation method for y such that

where F(x) is solely a function of x and G(t) is only a function of time. By substituting Eq. (S8.5) into Eq. (S8.4) (without the damping effects and magnetic actuation) and rearranging the terms, we obtain:

Based on classical vibration analysis, both sides of Eq. (S8.6) are independent of x and t, and they can be equated to become a function of the R ^th natural frequency, of the robot:

and

After solving the homogeneous equation in Eq. (S8.8), we can express the R ^th mode shape of the robot (F _R ) with ω _{η Κ}

where and the variables A _R , B _R , CR and D _R are constants dependent on the

robot's boundary conditions. Because the robot has free ends, we substitute the fol- lowing free-free boundary conditions into Eq. (S8.9):

Physically, the boundary conditions in Eq. (S8.10) dictate that there are no shear forces and bending moments acting on the free ends of the robots. Based on these boundary conditions, we can obtain the following characteristic equation for the robot:

Once the characteristic equation has been solved numerically, the robot's three lowest natural frequencies are given as: The first mode, ω _{η 1} , represents the robot's rigid-body motion, while the second and third mode, ah, 2 and oh,3, are the two lowest non-zero fundamental natural frequency of the robot. As the first mode does not describe the shape change of the robot with respect to the actuating magnetic field's frequency, we begin our frequency response analysis with ω _{η 2} and ω _{η 3} . Based on Eq. (S8.12), these two natural frequencies can be expresse

The mode shapes corresponding to ω _{η 2} and ω _{η 3} affect the small deflection cosine and sine shapes shown in Fig. 19b(ll-lll), respectively. Therefore, by using Eq. (S8.13), Eq. (S8.9) and Eq. (S8.5), the undulating traveling wave generated along the soft robot's body can be expressed as:

where the variables Ri and /¾ represent the asymptotic approximation of the Bode magnitude function of ideal 2 ^nd order systems, which can be expressed as:

To better analyze Eq. (S8.14), we rearrange its terms to become:

Physically, the first and second components on the RHS of Eq. (S8.16) represent the traveling and stationary waves along the robot's body, respectively. As the stationary wave component is a time-symmetrical motion, it does not propel the robot at low Reynolds number regime i.e., Re « 1 . Hence this component can be excluded, and the effective undulating swimming gait of the robot, i.e., the effective traveling wave, can be expressed as

Based on the Taylor's swimming sheet model, the swimming speed, Vswi _m , of this effective undulation swimming gait can be represented as

where f represents the traveling wave's frequency (coincident with the rotation frequency of B). The equation implies that: 1 ) when when

Because this locomotion can be represented as a low pass system, our scaling analysis will be based on the range of f smaller than ω _{η 2} .

Equation (S8.18) is valid for two types of boundary conditions: 1 ) the robot is swimming in the bulk of a fluid, and 2) the robot is swimming on an air-fluid interface. The key assumptions for the Taylor's swimming sheet model in Eq. (S8.18), however, is that it assumes low Reynolds number, i.e., Re « 1 , and also the wavelength is much larger than the traveling wave's amplitude, i.e.,

C. Discussion

To discuss how the undulating swimming speed, Vswim, varies against the robots' dimensions, we include the experimental data of six additional robots that have either a different L, w or h compared to the original robot (Fig. 19a). We conducted two types of experiments on these robots, to respectively evaluate their 1 ) V _SW im versus f, and 2) Vswim versus B. In all of these experiments, our robots were pinned at the air/water interface.

For the first type of experiments, we fixed B at 5 mT and varied f from 40 Hz to 160 Hz. The experimental results are shown in Fig. S26 together with the theoretical predictions based on Eq. (S8.18). The experimental data indicate that, for all robots, Vswim peaks at a robot-specific frequency, and these peak frequencies range between 60 Hz to 100 Hz. Our theoretical model predicts the existence of peak swimming speeds for the robots, though it is unable to capture the specific f at which the peak speed occurs. Despite this limitation, the model can approximate V _SW im versus f to be around one order of magnitude. The biggest discrepancy between the theoretical pre- dictions and experimental results is illustrated in Fig. S26c(ll), for which the experimentally obtained V _SW im is about 20 times faster than the predicted speed when f > 100 Hz.

For the second type of experiments, we varied B from 1 mT to 5 mT while fixing f at 40 Hz. The experimental results are shown in Fig. S27 together with the theoretical predictions based on Eq. (S8.18). Our experimental results show that V _SW im has a positive correlation with B and this agrees with our theoretical prediction. In general, the model can approximate V _SW im versus B within the same order of magnitude for all the robots.

The experimental data for the dependency of V _SW im on f and B parameterized by the dimensions of the robots are shown in Figs. S28 and S29. In general, these experimental data show that robots with larger L can swim faster and that w does not have an obvious correlation with V _SW im, and these phenomena agree with our theoretical predictions. However, in contrast to our theoretical predictions, our experiments indicate that robots with larger h will swim faster.

The discrepancies between the experimental results and theoretical predictions can be attributed to several factors. First, the Taylor swimming sheet model in Eq. (S8.18) assumes Re « 1 , whereas our characterization experiments displayed in Fig. S28 and S29 show that Re is to be between 4.6 and 190 for the original robot used in the main text, and between 2.54 and 676 for all the robots tested in SI. Since inertial effects are not included in the model, a discrepancy between the model and the data is expected. Indeed, to our best knowledge, precise solution of swimming at high Reyn- olds number usually should be solved numerically and it is non-trivial to include the inertia effects in our analytical model. Second, the soft robots violate the model assumption that the swimmer's wavelength is much bigger than its traveling wave's amplitude, as the amplitude and wavelength are in the same order of magnitude in some conditions. Third, uncertainties in our fabrication process cause our robots to have a residual small curvature even in its rest state (SI section 1 D and Fig. S1 e). Such uncertainty makes the amplitude of the sine shape smaller than predicted (see Fig. 19b(ll)). The ensuing non-ideal traveling wave across the robot body introduces an additional discrepancy between the model and the experimental results. Lastly, we did not consider surface tension effects in the model because the dynamic angles formed by water at the edges of the robots, which depend on the amplitude and frequency of the traveling wave as well as the material composition of the robots, are not trivial to be accounted for.

In summary, our experimental data suggest that robots with larger h and L can swim faster, while no obvious correlation between and w has been observed. For most cases, our theoretical model provides an order-of-magnitude accurate approximation for Vswim- Therefore, the model can be used as a guide for designers to estimate the swimming speed of future soft robots that can produce a similar swimming gait. We will try to develop a more accurate model that can predict at high Re regime and account for the robot's surface effects in the future. We will also further improve our fabrication method such that the robots can produce a better undulating traveling wave along their body.

S9 - Analysis for crawling

To analyze the crawling locomotion, here we specify the B sequence required to invoke this locomotion (section S9A), and based on our observations from the experi- ments conducted for seven soft robots with different dimensions (section S9B), we propose a suitable fitting model to account for it (section S9C).

A. Experimental details

The crawling locomotion shown in Fig. 20g was performed in a glass tunnel of rectangular cross-section (0.645 mm χ 2.55 mm) using a rotating B in the XY plane with frequency f of 20 Hz (counter-clockwise) and B of 30 mT (Fig. S30). The B sequence here can be described by Eq. (S8.1 ). While the crawling locomotion is possible for B > 2 mT for the original robot, we used B = 30 mT for the experiment in Fig. 20g to obtain a noticeable crawling displacement per cycle. Because the crawling robots were constrained by the walls, their deflections were small even when a large B was applied.

To better understand the crawling locomotion, we fabricated six additional robots, having either a different L, w or h compared to the original robot shown in Fig. 19a, and conducted two types of experiments to study respectively 1 ) crawling speed, Vcrawi, versus f (from 20 Hz to 140 Hz) while fixing B at 10 mT, and 2) Vcrawi versus B (from 2 mT to 18 mT) while fixing f at 20 Hz.

B. Discussion

Similar to the Taylor swimming gait presented in section S8, a rotating B is used for generating a traveling wave along the robots' body to enable the crawling locomo- tion (Fig. S30). The key difference between crawling and undulating swimming, however, is that the crawling direction is parallel to the propagation direction of the traveling waves while the motion of the swimming gait is anti-parallel to its traveling wave.

The experimental results in Fig. S31 show the variation of V _C rawi versus f parameterized by the dimensions of the robots. It suggests that V _C rawi of almost all robots in- creases linearly with f until f = 40 Hz. When f reaches 60 Hz, V _C rawi for the robot with h = 333 μιτι starts to decline. Such phenomenon affects also other robots when f becomes higher. We remark that this happens when the traveling wave can no longer produce constant amplitude: The robot is still able to crawl, but the relationship between Vcrawi and f is no longer linear.

The experimental data for the dependency of V _C rawi on B parameterized by the dimensions of the robots are shown in Fig. S32. In general, we observe a positive correlation between V _C rawi and B, until V _C rawi starts to saturate for relatively large B. Such saturation of V _C rawi could be explained by the change of the contact area between the robot and the tunnel as we vary B (Fig. S33). Using for instance the original robot (with dimensions of 3.7 mm χ 1 .5 mm χ 0.185 mm), we observe that V _C rawi starts to saturate after B exceeds 14 mT, and that this corresponds to a saturation of the contact length between the robot and the inner surface of the tunnel (Fig. S33). This is expected, since the amplitude of the traveling wave across the robot's body length within the tunnel is eventually constrained by the cross-section of the tunnel.

At last, we do not observe an obvious correlation between the robots' crawling speed and their dimensions L, w, or h: For any specific value of f and B, appears similar for all robots. C. Fitting model

The observations reported in the previous subsection suggest that Vcrawi can be approximated with a polynomial function of f and by an exponential function of B. We consequently provide the following fitting model:

where a, k ₀ , k _1t and k ₂ are the fitting parameters. Dimensions of the robots are not included in the model, since they did not show any obvious co-relationship with the experimental crawling speeds. To estimate the best fitting parameters, we formulated the following optimization problem:

where the subscript / ^' represents the corresponding ^h data point in the plots of Fig. S31 and Fig. S32. Numerical solutions of Eq. (S9.2) with MATLAB 2016 found the best fitting parameters to be a = 0.606, k ₂ = 0.005, / i = 0 and k ₀ = - 0.009. The Vcrawi predicted for each robot using the fitted model is plotted in Figs. S34 and S35. We remark that the applicability of this fitted crawling model is limited to soft robots that have dimensions of 185 μΐη≤ h < 333 μητι, 2.4 mm≤ L < 5 mm and 0.75 mm≤ w < 2 mm, crawling in a glass tunnel that has a rectangular cross-section of 0.645 mm χ 2.55 mm using a rotating B that has the following range of parameters: 2 mT≤ B < 18 mT and 20 Hz≤f≤ 140 Hz.

S10 - Analysis for jellyfish-like swimming

We discuss the jellyfish-like swimming by first specifying typical required se- quences of B to enable this locomotion (section S10A) and then describing our experimental results for seven robots of different dimensions (section S10B). Deriving a comprehensive analytical model for such complex swimming gait (Re > 1 ) is an open quest ¹⁴ and is beyond the scope of this paper. Therefore, here we provide a fitting model that can approximate the observed swimming speed of our robots (section S10C).

A. Experiment details

To perform the jellyfish-like swimming locomotion for the original robot shown in Fig. 20a, we used the sequences of B given in Fig. S36a in the XY plane, where the maximum B, B _max = 17 mT and f = 25 Hz. Theoretically, B should oscillate along a = 135° (equivalently to a = 315°) to create similar 'C'-shape and 'V'-shape configurations that are shown in Fig. 19b. Due to fabrication uncertainty, Fig. S36b however shows that the robot can perform this locomotion when B is oscillating along the direction a = 105° and a = 285°. We have therefore calibrated the applied B for this specific robot in our experiments. To characterize the jellyfish-like swimming, six additional robots with different L, w or h compared to the original robot were also fabricated. Like the jumping locomotion described in SI section S4, the characterization here is done with robots that have /3R = - 90°, instead of /3R = 45°as this allows us to simplify our setup to become a pair of Helmholtz coil. Using these experimental conditions, the robot can generate a jellyfishlike swimming locomotion when we prescribed a B sequence that has time-varying magnitude along the Y-axis of the global reference system. Mathematically, B can then be simply expressed as

Using these robots, we have conducted two types of experiments to evaluate re- spectively 1 ) the jellyfish-like swimming speed of the robots, Vj _f , versus f (from 20 Hz to 40 Hz) with Bmax fixed at 20 mT, and 2) j _f versus B _max (from 20 mT to 40 mT) with f fixed at 20 Hz.

B. Discussion

The prescribed B sequence (Fig. S36a) is meant to make the robot implement a slow recovery stroke where the robot changes from its 'C'-shape to 'V'-shape configuration (0 ms - 19.5 ms in Fig. 20a) and a fast power stroke (19.5 ms - 32 ms in Fig. 20a), which makes the robot change back to its 'C'-shape configuration again. Based on the Fig. S38 and S39, this swimming gait has Reynolds number that ranges from 74 to 190 for the original robot described in main text, and from 39 to 421 for all the robots tested in the SI. This gait resembles that of an actual jellyfish swimming, e.g., as presented in According to that study, part of the remarkable efficiency of jellyfish swimming is to be attributed to exploiting the vortex rings produced by the sequence of power and recovery strokes. We evidenced in Fig. S37 the generation of fluid vortices in the wake of our robots' jellyfish-like swimming locomotion by using 45 μιτι diameter polystyrene tracer microspheres (Polysciences, Inc).

Dependency of the experimental j _f on f is shown in Fig. S38, where the three sub-figures allow to better identify trends of the robots' swimming speed as we change their dimensions. As shown, except for the robot with L = 4.35 mm (Fig. S38b) and w = 2 mm (Fig. S37c), j _f peaks within 30 Hz < / < 40 Hz and starts to decay afterwards. Indeed, for/ > 40 Hz, the robots are observed to become unable to fully complete their recovery and power strokes within each period.

Dependency of the experimental j _f on B _max is shown in Fig. S39. In general, we observe that the robots swim faster as we increase B _max . However, j _f of more compliant robots start to decrease when B _max exceeds a robot-specific threshold value. For example, the robot with h = 108 μιτι starts to swim slower after B _max exceeds 25 mT (Fig. S39a). This arguably happens because once B _max > 25 mT the speed of the robot's recovery stroke becomes so fast that the robot is observed to produce significant downward motion. Such retarding effect is also observed in other robots, such as the robot with h = 142 μιτι in Fig. S39a and the robot with L = 5 mm in Fig. S39b. Based on the experimental data shown in Fig. S38 and S39, we observe that more compliant robots (i.e., larger L, smaller h) swim faster for B _max < 30 mT. This is expected, since these robots can create larger power strokes compared to their stiffer counterparts when subjected to the same B _max . Based on Fig. S38c and S39c, w only affects the dependency of j _f on f but not on B _max . Based on Fig. S38c, robots with larger w generally have a relatively lower j _f as we increase f.

C. Fitting model

The observations presented above suggest that we can fit polynomial functions for both f and B _max . Moreover, as robots with different dimensions have peak speeds for different values of f and B _max , the coefficients of the polynomial functions need be a function of h, L and w. In view of these observations, we propose the following fitting model:

where k a,, b and c, are the fitting parameters. To estimate the best fitting parameters, we formulate the following optimization problem:

where t e su scr pt represents t e correspon ng ata po nt n t e p ots presented in Fig. S38 and Fig. S39. By solving Eq. (S10.2) numerically with MATLAB 2016, we present the obtained fitting parameters in Table S3. The fitting model in Eq. (S10.1 ) shows that the relationship between Vj _f and the robot's dimensions is highly nonlinear.

The values of j _f predicted for each robot by the fitted model are shown in Fig. S40 and Fig. S41 . We remark that the applicability of this fitted swimming model is for soft robots that have dimensions of 108 μηη≤ h < 185 μητι, 3.7 mm≤ L < 5 mm and 1 .5 mm ≤ w < 2.5 mm, using a sequence of B that has the following range of parameters: 20 mT≤ B _m ax≤ 40 mT and 20 Hz≤ f ≤ 40 Hz.

S11 - Transition modes

A. Immersion

Driven by the B sequence shown in Fig. S42a in the XY plane, the immersion shown in Fig. 20d is realized by the following steps:

1 . The initial B (B = 20 mT, a = 90°) is imposed to curl the robot downward (32 ms in Fig. 2d). 2. The robot rotates counterclockwise to reduce body contact with the water/air interface to a minimum (41 ms in Fig. 20d), leaving a single end of the robot pinned to the water surface.

3. A quick 180° flipping of the external B disengages the pinned end of the robot from the water surface (66 ms in Fig. 20d, corresponding to 0.86 time unit in Fig. S42).

Consequently, the robot, denser than water, sinks.

We remark that for our robots, with main dimension L « L _c , the Bond ber therefore the magnitude of surface tension effects cannot be neglected compared to gravitational ones. Immersion in the liquid bulk re- quires that the robots break the surface tension of water. The minimum force F _im necessary to break the contact between a single robot end and the water/air interface (as for the robot configuration in Fig. 3b) is of the order of γρ, where p = 2(w + h) is the perimeter of the cross-section of the robot and γ = 72 mN-m ^"1 is the surface tension of pure water at room temperature. In the case of the original robot, F _im « 243 D243 he surface tension of pure water at room temperature. In tB in step 3 above.

B. Landing

Driven by the B sequence shown in Fig. S43 in XY plane, the transition from the top of the water meniscus to the adjacent platform shown in Fig. 20c is realized by the following steps:

1 . The initial B {B =1 1 .7 mT, a = 270°) is imposed to curl the robot (2.8 s in Fig. 20c, 0 time unit in Fig. S43a). During this motion, water can easily dewet the robot surface (see next subsection).

2. Subsequently the B starts to rotate clockwise. The B magnitude decreases to 4.7 mT at the 0.26 time unit in Fig. S43a since most of the robot has peeled away from the water surface and a large B is no longer necessary.

Here we want to highlight that, due to residual strain energy caused by fabrication uncertainty (section S1 D), this specific robot in Fig. 20c can perform rigid body rotation even when B ceases to be 4.7 mT. We have therefore calibrated the applied B for this specific robot in our experiments.

C. Effect of robot's surface properties

The native elastomer surface of the robot is hydrophobic and microscopically rough (SI section S1 B). The low wettability of the robot surface, and particularly the high receding water contact angle (-78°, see SI section S1 B), leaves the receding triple contact lines free to move during the rotation of the robot. Hence, the robot can easily peel away from the water surface to stand on the adjacent solid substrate (Fig. 20c). Surface hydrophobicity also enhances the pinning of the robot at the water/air interface, because it extends the angular range that can be spanned by the edge angles without displacing the position of the contact line on the robot body. In combination with its microscale roughness, surface hydrophobicity promotes the capture and adhe- sion of microscopic air bubbles on the robot surface upon immersion in water from air. The presence of a plastron, i.e., a sparse layer of microbubbles, is revealed by the shiny appearance of the robot underwater, caused by the light refracted at the water/air interface of the bubbles. The plastron marginally increases the water volume displaced by the submerged robot, and thus its buoyancy.

Moreover, surface hydrophobicity energetically promotes the emersion of the robot from the water bulk to the water surface, since the replacement of the elastomer/water interface with a full elastomer/air interface during emersion decreases the interfacial energy of the system. On the other hand, for the same reason, hydrophobicity also makes it harder for the robot to disengage from the water surface during immersion. S12 - Scaling analyses summary

A. Scaling analysis for robot's physical parameters

Here, we provide a simple scaling analysis for the robot's physical parameters. First, the magnitude of the magnetization profile for the infinitesimal elements shown in Fig. S3 scales proportionally with h and w, while the robot's mass and the magnitude of net magnetic moment, M _net , scales proportionally with h, L and w. Second, the second moment of inertia of the robot, /, has the following scale laws: / oc h ³ and / oc w. The parameter / will be important for describing the bending stiffness of the robot.

B. Scaling analyses for locomotion modalities

Based on the modeling and experimental results discussed in SI sections S4-10, this sub-section and Table S4 will provide an overview of scaling analyses for each locomotion modality. We also suggest the operational range of f for each locomotion modality. As long as f is within this range, all of our tested robots are able to execute their desired locomotion successfully. Table S4 also includes the best performance index observed from each locomotion modality.

We remind that our analyses and discussions are intended to facilitate future studies on soft-bodied locomotion. In the future, we will also expand these analyses such that they can be used for optimizing the performance of miniature robots with multi- locomotive capabilities.

1 . Jumping. Performance index: maximum jumping height, H _max .

Theoretically, H _max is independent of w, and its relationships with L and h are highly non-linear. The theory, however, still suggests that H _max has a positive correlation with L and a negative correlation with h. Experimentally, we observe that robots with larger L and thinner h can jump higher. The data do not show observable trends for w. Hence, these experimental results agree with the trend predicted by the theory.

Rolling. Performance index: rolling speed, V _ro n.

Theoretically, ν _Γθ11 is independent of w and h when the rolling robots have similar curvatures. By using simple approximations, we can deduce that V _roll oc L. Experimentally, good control over rolling is achieved for all the robots when f is between 0 Hz and 4 Hz (Fig. S12). Our data also suggest that when the robots have similar curvatures, robots with larger L will roll faster. There are no observable effects for w and h. These experimental results agree with the trend predicted by the theory.

Walking. Performance index: walking speed, V _wa ik-

Theoretically, V _walk is independent of w and its relationship with L and h is nonlinear. However, the theory does imply that V _walk has a positive correlation with L and a negative correlation with h.

Experimentally, good control over walking is achieved for 2 Hz < ^" < 1 1 Hz (Figs. S17). Robots with larger L and thinner h are observed to walk faster when subjected to the same S _max (Fig. S18). The data do not show observable trends for w. These experimental results agree with the trend predicted by the theory.

Meniscus climbing: Performance index: minimum required B, i.e., B _mm to reach the top of the meniscus.

Theoretically, 5 _min oc L ^~3 , h ³ . The current model does not account for the effects of w.

Experimentally, all the tested robots can climb to the top of the water meniscus shown in Fig. S20 for S _min > 6 mT (Fig. S23). Robots with a larger L and thinner h require a smaller S _min to ascend to the top. These experimental results agree with the trend predicted by the theory.

Undulating swimming: Performance index: swimming speed, V _swim .

Theoretically, we have: 1 ) when when

As the robotic system is a low-pass system, we will consider only

the first condition (i.e., when for the scaling analysis. Vswim is always

independent of w.

Experimentally, all the tested robots can swim successfully when f is between 40 Hz and 160 Hz (Fig. S28). Robots with larger L and thicker h can swim faster (Fig. S29). The experimental data do not show observable trends for w. There exists a discrepancy between the experimental and theoretical scaling for h.

6. Crawling: Performance index: crawling speed Vcrawi-

Experimentally, good control over crawling is achieved for 20 Hz < / < 40 Hz (Figs. S31 ). The data suggests that V _C rawi has no observable dependency on L, w, and h.

Based on the experiment data, the fitting model suggests that V _C rawi is independent of w, L, and h.

7. Jellyfish-like swimming: Performance index: swimming speed Vjf.

Experimentally, all the tested robots can swim successfully for 30 Hz < / < 40 Hz (Figs. S38). It is observed that Vjf may peak at certain B _max (Fig. S39), which depends on specific robot design. After exceeding this B _max , the robots may produce very fast recovery strokes that reduce Vjf. Before reaching this B _max , we observe that robots with larger L and thinner h swim faster. The dimension w only affects the dependency of Vjf on f, but not on B: robots with larger w generally have a relatively lower Vjf as we increase f.

Based on the experimental data, the fitting model suggests that the Vj _f has a highly non-linear relationship with L, h, and w. The values for the fitting parameters are provided in Table S3.

In the dimension range we have investigated, our theoretical and fitting models predict that a larger L and a smaller h are always preferred for multimodal locomotion as it can help the robot to move faster and jump higher. The models also suggest that w will only affect the jellyfish-like swimming locomotion and minimizing w can help to increase this swimming speed.

C. Upper and lower bounds for L, w, and h

Based on our models discussed above, the robot may not be able to fulfill all modes of locomotion if its L is too small, or when its h and w are too large. In addition, the upper bound of L, w, and h is typically constrained by the size requirements of specific applications as well as the maximum allowable workspace of the electromagnetic coil setup that generates the spatially uniform B. The lower bound of h and w is set by our current fabrication limits where we could only demold robots with h > 40 μηη and w > 0.3 mm repeatedly. The lower bound of L is 1 mm because it is very challenging to manually wrap sub-millimeter beams onto the circular jigs, shown in Fig. S1 a, during the magnetization process. S13 - Multimodal locomotion

Here, we describe the procedures to implement the robot's multimodal locomotion shown in Fig. 3. In Fig. 3a-d, the environment is composed of laser-milled, poly(methyl methacrylate) platforms. The platform is treated with oxygen-plasma and has a water contact angle < 10°. The platform is encased in a Plexiglas box (52 mm χ 32 mm χ 25 mm). The liquid in this environment is de-ionized water.

In Fig. 3a, the robot first uses the rolling locomotion to fall onto the water surface. After dipping into water, the robot drifts slightly away from the platform along the meniscus and starts swimming after it has stabilized. A clockwise B (B = 3 mT, 25 Hz) is applied to let the robot swim rightwards. In Fig. 3b, the robot sequentially bends downwards on the water surface and rotates counterclockwise to sink into the water, and subsequently it uses the jellyfish-like locomotion to swim back to the water surface in Fig. 3c.

Figure 3d shows the robot using four locomotion modes to navigate through liquid and solid surfaces. After climbing onto the solid substrate, the robot uses the directional jumping strategy to overcome the otherwise unsurmountable obstacle. By exploiting the robot's net magnetic moment, we rotate the robot to its required initial pose for implementing the directional jumping. Once the robot has jumped across this obstacle, it starts to walk. The robot 10 incurs slight slipping before and after jumping across the obstacle. This slipping is induced by unwanted magnetic gradient-based pulling forces generated by the spatial gradients of B, which are present because the robot has exceeded the 95% homogeneous region that can provide spatially-uniform B for the robot (see SI section S2). We also remark that this slipping generally happens when the robot has only one end in contact with the substrate. Since the friction is the smallest at this moment, even a small magnetic gradient-based pulling force can translate the robot. However, while there exists such an external disturbance for our robot, it is worth noting that the robot can still walk successfully even against the gradient after jumping over the obstacle. Moreover, since the deformed robot has a net magnetic moment, we are able to steer the robot along a desired direction even after it was disrupted by the magnetic gradient-based pulling forces (see SI section S3B).

In Fig. 3e, the robot first walks towards a glass tunnel (inner diameter: 1 .62 mm), which eventually blocks its path. As the tunnel is too small for the robot to walk in it, the robot switches to its crawling locomotion by changing the magnetic field input to a rotating B. After the robot has crawled through the tunnel, it resumes its walking loco- motion until it leaves the working space. S14 - Toward biomedical applications

A. Locomotion in a surgical phantom

Here we demonstrate that our robot can use a combination of meniscus climbing, landing, rolling and jumping to fully explore a surgical human stomach phantom (Fig. 21 a).

The robot moves back to the starting point very quickly, because it is pulled by unwanted magnetic gradient-based pulling forces generated by the spatial gradients of

B, which are present because the robot is outside the region that can provide spatially- uniform B (see SI section S2). While there exists such an external disturbance for our robot, it is worth noting that the robot can still roll successfully under such conditions.

B. Ultrasound-guided locomotion

Here we demonstrate that our robot has the potential to be integrated with biomedical imaging systems towards realizing potential in-situ applications like minimally invasive surgery. In particular, we demonstrate ultrasound-guided locomotion (Fig. 21 b, Fig. S44). For the experiments, we have used three different biological phantoms to mimic biomedical landscapes, i.e., a water reservoir, tunnels and gaps. We inserted the ultrasonic scanning head of the system (Fujifilm, Vevo3100 Imaging System) from the top of our magnetic setup and put it in contact with the phantoms through a thin layer of ultrasonic gel. The ultrasound system was found able to track the position of our robots while they were performing jellyfish-like swimming, crawling, and rolling within the narrow spaces concealed by chicken muscle tissue. This allowed us to image and track the robots as they navigated in the phantoms. We were able to produce clear images of the robots while having no disruptions to our actuation system.

We also report that while the technicians who worked with us on the ultrasonic setup can easily detect the motions of the robot from the ultrasound images, it is possible to develop imaging algorithms to automate this process in the future.

C. Cargo delivery

The design of the modified robot, presented in Fig. 21d, for selective, magnetically-triggered drug release, presented is shown in Fig. 18. An extra strap is added to the original robot body. For clarity, we only mark out the m _y component of the strap's magnetization profile, as it is the dominant and functional component in releasing the drug. During locomotion, the cargo is mechanically bound on the robot by inserting the strap head into a hole on the robot's body. When a large B along the y-axis, B _y (-19 mT), is applied, the magnetic torque bends and unlocks the strap to release the cargo. Due to the restriction of the hole, a smaller S _y (<18 mT in our case) is not strong enough to open the strap and hence this could allow us to preserve our locomotion modes. S15 - Spatial gradient of B

By creating a spatially non-uniform B, we can exert magnetic forces on miniature devices that have a net magnetic moment M _ne t. The mathematical description for these applied forces can be expressed as:

The method that uses the spatial gradients of B to translate the miniature robots is typically known as magnetic gradient-based pulling or gradient pulling ^{1 ,5} . As our robots can possess an effective M _ne t when they deform, we can also use the spatial gradients of B to exert magnetic gradient-based pulling on the robots such that their locomotion capabilities can be enhanced. For instance, we can use magnetic gradient-based pull- ing to modulate the robot's jump height for the jumping locomotion.

Theoretically, magnetic gradient-based pulling may also levitate a magnetic robot and make it traverse across different terrains and obstacles. However, this approach is not practical as the dynamics of this actuation mode is inherently unstable, especially if the small-scale robot has to be levitated against gravity. Furthermore, Abbot et al. has shown that time-asymmetrical locomotion like the undulating swimming gait is more energy efficient than simple magnetic gradient pulling actuation methods, i.e., they produce higher propulsion forces. It is also easier to actuate miniature magnetic robots using B than its spatial gradients when their distance from the external electromagnets increases. As many biomedical applications require the external electromagnets to be placed outside of the human body, the considerable distance between the electromagnets and the proposed robots would make magnetic gradient-based pulling methods less appealing.

Supporting Figures also appended to the application

Fig. S1 . Magnetization process of the soft magneto-elastic robot, the experimental results for the magnetization magnitude and Young's modulus of similar magnetic elastomers with differing mass ratios, and the experiment for showing different rest state residual curvatures, (a) Sketch of the elastomeric beam, loaded with NdFeB mi- croparticles, wrapped around a cylindrical glass rod (perimeter: 3.7 mm) and magnetized by a uniform 1 .65 T magnetic field B produced by a vibrating sample magnetometer (VSM). The relative orientation _R determines the phase shift in the resulting magnetization profile, (b) The programmed magnetization profile in the unwrapped robot body in its local xyz frame ( R = 45°). (c) Experimental values of the magnetization magnitude for soft magnetic materials with different mass ratios, (d) Experimental values of E for soft magnetic materials with different mass ratios. The mass ratio is defined as the ratio of mass of NdFeB particles to the mass of Ecoflex-10, e.g., 0.25:1 , 0.5:1 , and we normalize these ratios such that it is always one part on the Ecoflex-10 side. Each data point in (c) and (d) was evaluated with three samples (n = 3). (e) The robot has two rest state residual curvatures (shown in III, V), which can be alternated via shape change. The five shapes of the robot are achieved sequentially by subjecting the robot to the B field sequence specified below the snapshots.

Fig. S3. Quasi-static analysis of the soft robot. Sketch of the geometry is not to scale for representation purposes. The bending moment acting on an infinitesimal element of the robot under steady-state deformation is shown in the inset.

Fig. S4. Soft robots implementing jumping locomotion, (a) Directional jumping by a soft robot with _R = 45° (also shown in Fig. 20h). (b) Straight jumping by a soft robot with = -90°. The magnetization profiles of the robots shown in (a,b) and their corresponding responses to the applied field are sketched in (c). Scale bars: 1 mm.

Fig. S5. The required B for the directional jumping locomotion, (a) Time sequence of B. B _x and B _y represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) shows the initial shape of the soft robot before it jumps (D R = 45° and reproduced from Fig. 20h). Scale bar: 1 mm.

Fig. S6. Sketch for the analysis of the straight jumping locomotion (not to scale). Fig. S7. Experimental and model-predicted results for H _max in the straight jumping locomotion, (a) Results for the original robot (dimensions defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation of three trials. Error bars in model-predicted values reflect the dispersion in image analysis (n = 3).

Fig. S8. The compiled experimental data, displaying how the robot's dimensions will affect Hmax for the straight jumping locomotion (H _max versus B). (a) Complied experimental results for robots that have different h. (b) Complied experimental results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials.

Fig. S9. The required B for the rolling locomotion, (a) Time sequence of the rotating B. B _x and B _y represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshot shows the initial 'C-shape' configuration of the robot before it rolls (reproduced from Fig. 20e). Scale bar: 1 mm.

Fig. S10. Sketch for the analysis of the rolling locomotion (not to scale). The illustration shows the robot (a) before and (b) after it has rotated a rotational displacement of §R. The center of mass is marked out by the black dot. This sketch also shows the applied forces and torques on the robot as it rolls.

Fig. S1 1 . Experimental and model-predicted results for the rolling locomotion (V _ro n ver- sus f). The applied B for all the robots is adjusted in such a way that all these robots can produce a curvature similar to that of the original robot when it is subjected to B =18.5 mT. (a) Results for the original robot (defined in Fig. 19a), where the applied B = 18.5 mT. (b) Results for robots that have a different h from the original robot. (I) B = 6.3 mT (II) B = 37 mT. (c) Results for robots that have a different L from the original robot. (I) B = 44 mT (II) B = 10 mT. (d) Results for robots that have a different w from the original robot. B = 18.5 mT for both (I) and (II). Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S12. The compiled experimental data, displaying how the robot's dimensions will affect Vroii for the rolling locomotion (V _ro n versus f). B is adjusted to make all the soft robots have similar curvature as that of the original one when it is subjected to B =18.5 mT. (a) Compiled experimental results for robots that have different h. B = 18.5 mT (h =185 μηη), B = 6.3 mT (h =108 μηη). B = 37 mT (h = 264 μηη). (b) Compiled experimental results for robots that have different L. B = 18.5 mT (L = 3.7 mm), B = 44 mT (L = 2.4 mm), B = 10 mT (L = 5 mm), (c) Compiled experimental results for robots that have different w. B = 18.5 mT for all the three cases. Each experimental data point presents the mean and standard deviation out of three trials.

Fig. S13. The required B for the walking locomotion, (a) Time sequence of B. This sequence can be sub-divided into three phases, i.e., 1 ) tilt forward, 2) tilt backwards, and 3) extend its front end. B _x and B _y represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshot shows the 'C- shape configuration during it walks (reproduced from Fig. 20f). Scale bar: 1 mm.

Fig. S14. (a) Stretched and (b) curled state of the robot during a walking cycle (not to scale). Fig. S15. Experimental and model-predicted results for the walking locomotion (V _wa ik versus f while B _max is fixed at 10 mT). (a) Results for the original robot (defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S16. Experimental and model-predicted results for the walking locomotion (V _wa ik versus B _max while f is fixed at 5 Hz), (a) Results for the original robot (defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Experimental data points present mean and standard deviation out of three trials (n = 3). Fig. S17. The compiled experimental data, displaying how the robot's dimensions will affect Vwaik for the walking locomotion (V _wa ik versus f while B _max is fixed at 10 mT). (a) Compiled experimental results for robots that have different h. (b) Compiled experimental results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S18. The compiled experimental data, displaying how the robot's dimensions will affect Vwaik for the walking locomotion (V _wa ik versus B _max while f is fixed at 5 Hz), (a) Compiled experimental results for robots that have different h. (b) Compiled experi- mental results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S19. The required B for the meniscus climbing locomotion, (a) Time sequence of B. B _x and B _y represent the magnitude of the B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshot shows that the robot in a 'C'-shape configuration (for B with□ = 315°) displaces water from the water-air interface, as evidenced by the air engulfed within its body, while remaining pinned to the water surface. Scale bar: 1 mm.

Fig. S20. Meniscus climbing by a magnetic soft robot that has an L of 5 mm. The sequential snapshots show side-views of selected, stable positions assumed by the robot along a positive water meniscus for increasing values of body curvature induced by the external B. The robot is fully immersed in water, its shape conforms to the local interface profile, and the perimeter of its top surface is pinned at the water-air interface. Fig. S21 . Sketches of the geometries for the meniscus-climbing models, a) Geometry for the arc-shaped meniscus-climbing model. CG and CB represent respectively the center of gravity of the immersed robot and the center of buoyancy of the displaced water volume, b) Geometry for the cosine-shaped meniscus-climbing model. The model specifically describes the soft robot in close proximity of the adjacent wall at the top of the positive water meniscus, and assumes coincidence of gravity and buoyancy centers. The sketches are not to scale for illustration purposes.

Fig. S22. Results for the meniscus-climbing locomotion (robot curvature versus robot dimensions). Comparison of experimental and model-predicted values of curvature radius for soft robots climbing the entire height of the water meniscus (i.e., for d→ 0, Fig. S21 a).

Fig. S23. Results for the meniscus-climbing locomotion (B _min as a function of robot dimensions). Comparison of experimental and model-predicted values of B _min required for soft robots to climb the entire height of the water meniscus (i.e., for d→ 0, Fig. S21 a).

Fig. S24. The required B for the undulating swimming locomotion, (a) Time sequence of B (f = 40 Hz, B = 3 mT). B _x and B _y represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshot shows the cor- responding small-deformation shapes of the original robot (Fig. 19b). Scale bar: 1 mm.

Fig. S25. Vibrational analysis of a soft robot. Sketch of the geometry highlighting forces and torques acting on an infinitesimal element of the robot (not to scale).

Fig. S26. Experimental and model-predicted results for the undulating swimming locomotion (Vswim versus f while B is fixed at 5 mT). (a) Results for the original robot (de- fined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S27. Experimental and model-predicted results for the undulating swimming lo- comotion (V _swim as a function of B while f is fixed at 40 Hz), (a) Results for the original robot (defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S28. The compiled experimental data, displaying how the robot's dimensions will affect Vswim for the undulating swimming locomotion (V _swim as a function of f while B is fixed at 5 mT). (a) Compiled experimental results for robots that have different h. (b) Results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S29. The compiled experimental data, displaying how the robot's dimensions will affect Vswim for the undulating swimming locomotion (V _swim as a function of B while f is fixed at 40 Hz), (a) Compiled experimental results for robots that have different h. (b) Compiled experimental results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S30. The required B for the crawling locomotion, (a) Time sequence of B. B _x and By represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) shows the instants of the robot when it creates the deformed sine (top) and the cosine (bottom) shapes during a crawling cycle. Scale bar: 1 mm.

Fig. S31 . The compiled experimental data, displaying how the robot's dimensions will affect Vcrawi for the crawling locomotion (V _cra wi as a function of f while B is fixed at 10 mT). (a) Compiled experimental results for robots that have different h. (b) Compiled experimental results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials (n = 3). Fig. S32. The compiled experimental data, displaying how the robot's dimensions will affect Vcrawi for the crawling locomotion (V _cra wi as a function of B while f is fixed at 20Hz). (a) Compiled experimental results for robots that have different h. (b) Compiled experimental results for robots that have different L. (c) Compiled experimental results for robots that have different w. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S33. Crawling speed and contact length between the original robot and the inner tunnel surface as a function of B (f is fixed at 20 Hz). The contact length is defined in the top subfigure, and is assumed uniform across the width of the robot, i.e., contact area = contact length χ w. The flat lining of the contact length correlates with the saturation of the crawling speed. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S34. Experimental and model-predicted results for the crawling locomotion (V _cra wi as a function of f while B is fixed at 10 mT). (a) Results for the original robot (defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S35. Experimental and model-predicted results for the crawling locomotion (V _cra wi as a function of B while f is fixed at 20 Hz), (a) Results for the original robot (defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S36. The required B for the jellyfish-like swimming locomotion, (a) Time sequence of B. B _x and B _y represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshots show the 'C- and 'V'-shape configurations assumed by the robot before the recovery stroke (top) and before the power stroke (bottom), respectively. Scale bar: 1 mm.

Fig. S37. Fluid vortices produced by the power stroke during the jellyfish-like swimming locomotion. The video snapshots in (a) and (b) respectively show the front and lateral views of two experiment trails recorded separately. The flow was traced using polystyrene particles (diameter: 45 μιτι). The original robot (dimensions defined in Fig. 19a) was used for this experiment. 40 frames recorded at 2000 fps (corresponding to a time period of 20 ms) were overlapped to produce the images. Scale bars: 1 mm.

Fig. S38. The compiled experimental data, displaying how the robot's dimensions will affect Vj _f for the jellyfish-like swimming locomotion (V _jf as a function of f while B _max is fixed at 20 mT). (a) Compiled experimental results for robots that have a different h from the original robot, (b) Compiled experimental results for robots that have a different L from the original robot, (c) Compiled experimental results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3). Fig. S39. The compiled experimental data, displaying how the robot's dimensions will affect Vj _f for the jellyfish-like swimming locomotion (V _jf as a function of B _max while f is fixed at 20Hz). (a) Compiled experimental results for robots that have a different h from the original robot, (b) Compiled experimental results for robots that have a different L from the original robot, (c) Compiled experimental results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S40. Experimental and model-predicted results for the jellyfish-like swimming locomotion (V _j f as a function of f while B _max is fixed at 20 mT). (a) Results for the original robot (dimensions defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each exper- imental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S41 . Experimental and model predicted results for the jellyfish-like swimming locomotion (V _j f as a function of B _max while f is fixed at 20 Hz), (a) Results for the original robot (dimensions defined in Fig. 19a). (b) Results for robots that have a different h from the original robot, (c) Results for robots that have a different L from the original robot, (d) Results for robots that have a different w from the original robot. Each experimental data point presents the mean and standard deviation out of three trials (n = 3).

Fig. S42. The required B for the immersion locomotion, (a) Time sequence of B. B _x and By represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshot shows the initial configuration of the robot for the immersion (reproduced from Fig. 20d). Scale bar: 1 mm.

Fig. S43. The required B for the landing locomotion, (a) Time sequence of B. B _x and B _y represent the magnitude of B projected along the X- and Y-axis of the global coordinate system, (b) The video snapshot shows the initial configuration of the robot for the landing (reproduced from Fig. 20c).

Fig. S44. Ultrasound-guided locomotion. Here, the dimensions of the robots are similar to the original robot (Fig. 19a), i.e., L = 3.7 mm, w = 1 .5 mm and h = 0.185 mm. The robot is marked out by red dotted line in a-c. (a) The robot with _R = - 90° uses the jellyfish-like swimming locomotion with the following parameters of B sequence: f = 30 Hz and B _max = 30 mT. This experiment is conducted in a phantom where the bottom and top parts of the swimming reservoir were made by chicken muscle tissue, (b) The robot with _R = 45° crawls by using the following parameters of rotating B sequence: f = 10 Hz and B = 10 mT. This experiment was performed in a silicone tube (inner diameter of 2 mm) buried inside the chicken muscle tissue (see (e) for an enlarged im- age), (c) The robot with _R = 45° rolls with the following parameters of rotating B sequence: f = 2 Hz and B = 15 mT. This is a re-plot of Fig. 21 b. The rolling experiment is conducted in the gap created between two layers of chicken muscle tissue, (d) The ultrasound scanning head and the chicken muscle tissue phantom when they are placed in the electromagnetic setup shown in Fig. S2. (e) An enlarged image of the robot when it is within the buried tube. Scale bars are: 1 mm (a-c) and 1 cm (d).

Supporting Tables

Table S1. Nomenclature.

Table S2. Physical properties of the robot.

Table S3. Summary of the fitting parameters for the jellyfish-like swimming locomotion.

Table S4. Scaling analyses summary. We provide the performance indicator, best observable performance index, operating range of f, and robot dimension effects for each locomotion mode. The suggested operating range off is based on our tested robots.