3D-printed low-voltage-driven ciliary hydrogel microactuators

Single-layer microelectrodes fabrication

Extended Data Fig. 2 illustrates the fabrication process of single-layer microelectrodes. Below are the step-by-step details:

Step 1. Preparation of the polyimide substrate (Extended Data Fig. 2a). This work uses PI2611 (HD MicroSystems) as the backbone material. PI2611 is poured onto a blank glass substrate and then spin-coated at 1,500 rpm for 30 s. The glass substrate used to hold the printed hydrogel microactuators must have a thickness of less than 300 µm, as this is the working range of the TPP laser. For this purpose, we use a 180-µm-thick glass substrate for hydrogel printing. No specific thickness is required for the glass substrate used for moulded hydrogel fabrication.

Step 2. Curing the PI2611 polyimide substrate (Extended Data Fig. 2b). We place the spin-coated substrate on a hotplate and heat it from room temperature to 150 °C. We hold this temperature for 10 min and then increase it to 200 °C. The temperature ramp rate is 20 °C per minute. This temperature is maintained for 5 h to fully cure the polyimide.

Step 3. Photoresist coating (Extended Data Fig. 2c). We pour the positive photoresist AZ ECI 3012 (MicroChemicals GmbH) onto the polyimide substrate and spin-coat the photoresist for 30 s at 5,000 rpm.

Step 4. Soft baking (Extended Data Fig. 2d). We bake the positive photoresist at 90 °C for 90 s.

Step 5. Ultraviolet exposure (Extended Data Fig. 2e). We expose the substrate for 8 s using the MJB4 mask aligner (SUSS MicroTec). The ultraviolet density of this machine is 14.3 mJ cm−2 and the patterns are defined using a photomask.

Step 6. Post-exposure bake (Extended Data Fig. 2f). We bake the exposed substrate at 110 °C for 90 s.

Step 7. Development (Extended Data Fig. 2g). We develop the substrate in AZ 726 (MicroChemicals GmbH) developer for 60 s to reveal the patterns.

Step 8. Platinum sputtering (Extended Data Fig. 2h). We deposit a 150-nm-thick platinum (Pt) layer onto the substrate using sputtering.

Step 9. Lift-off process (Extended Data Fig. 2i). We dip the substrate in TechniStrip Micro D350 (MicroChemicals GmbH) photoresist stripper to remove unwanted material and obtain the final microelectrode structures.

Hydrogel solutions for 3D printing

Four hydrogel precursor solutions are prepared for TPP-based 3D printing. Acrylic acid (AAc, Sigma-Aldrich) and acrylamide (AAm, Sigma-Aldrich) served as monomers, N,N′-methylenebisacrylamide (BIS, Sigma-Aldrich) as the crosslinker and Omnirad 2100 (IGM Resins) as the photoinitiator. Ethylene glycol was used as the base solvent.

The four solutions differ in the mass fractions of AAc, AAm and the photoinitiator. The molar masses of AAc (72.06 g mol−1) and AAm (71.08 g mol−1) are nearly identical; therefore, maintaining a constant total mass fraction of AAc and AAm effectively corresponds to a constant overall molar concentration of monomers. This ensures that, across all formulations, the total monomer concentration remains comparable when only the relative ratio of AAc to AAm varies.

Because the printability of the precursor solution depends on the AAc-to-AAm ratio, the photoinitiator fraction is experimentally adjusted to achieve consistent printing quality. For each monomer ratio, several photoinitiator concentrations are tested under identical printing parameters and the optimal concentration is determined as the one that produced structures without overexposure or underexposure. Overexposure and underexposure are defined as follows: when printing a 10-µm feature under identical parameters, a fabricated structure larger than 10 µm is considered overexposed, whereas one smaller than 10 µm is considered underexposed. All mass fractions reported below are calculated with respect to the initial mass of ethylene glycol before the addition of monomers, crosslinker or photoinitiator.

Solution 1: AAc 15 wt%, AAm 60 wt%, BIS 2 wt%, photoinitiator 7.5 wt%.

Solution 2: AAc 30 wt%, AAm 45 wt%, BIS 2 wt%, photoinitiator 18 wt%.

Solution 3: AAc 45 wt%, AAm 30 wt%, BIS 2 wt%, photoinitiator 28 wt%.

Solution 4: AAc 60 wt%, AAm 15 wt%, BIS 2 wt%, photoinitiator 38 wt%.

Hydrogel microactuator 3D printing

A commercial TPP system (Photonic Professional GT, Nanoscribe GmbH) is used to fabricate hydrogel microactuators. The printing is performed using a 25× objective lens with an exposure power of 15 mW. The slicing and hatching distances are set to 300 nm and 200 nm, respectively, with a 45° hatching angle between adjacent layers (Extended Data Fig. 3).

We immerse the printed sample in an ethylene glycol bath to develop the hydrogel microactuator structure. Subsequently, the sample is transferred to a DI water bath for 10 min. We repeat this process three times to fully replace the solvent inside the hydrogel from ethylene glycol to DI water. Afterwards, the aqueous environment can be adjusted for different experiments accordingly.

Comparison between the centimetre–millimetre hydrogel and the micrometre hydrogel

In contrast to the fast bending mechanism of the micrometre hydrogel actuator, the internal ion migration is much slower in previously reported centimetre–millimetre-scale ionic hydrogels not fabricated by TPP49,50. Even at the same electric-field intensity as our gel microcilia system, for example, 10,000 V m−1, the H+ ions and Na+ ions would take 0.3 s and 16.7 s, respectively, to traverse a 1-cm distance, which is a typical thickness of a centimetre-scale hydrogel actuator. This ion migration time is orders of magnitude slower than at the micrometre scale. Moreover, its larger pore size (Extended Data Fig. 1b (i)) and smaller effective EDL surface area reduce ion transport, thereby weakening bending. Finally, achieving such high electric fields in large systems needs impractical voltages (for example, 200 V for 2-cm-spaced electrodes), which can trigger electrolysis and other vigorous electrochemical reactions.

In fact, the internal ion migration in millimetre-scale hydrogel is trivial compared with ion partitioning at the gel–liquid interface. Therefore, previously reported millimetre-scale hydrogels49,50 operate through mechanisms distinct from those proposed in the present work. They bend by means of bath-induced pH gradients or interfacial osmotic effects, driven by electrolysis or ion partitioning across the gel–solution interface42,43,44. By contrast, micrometre-scale hydrogels depend on internal ion migration and nanometre-scale channels, resulting in opposite bending directions and response times more than two orders of magnitude faster (Extended Data Fig. 5 and Supplementary Videos 6 and 7). For example, in DI water, a millimetre-scale actuator bends towards the anode in about 120 s, whereas a micrometre-scale actuator bends towards the cathode in about 0.2 s (Extended Data Fig. 5a,b). In 0.15380 mol l−1 NaCl, the millimetre-scale actuator bends towards the cathode in about 30 s, whereas the micrometre-scale actuator bends towards the anode in about 0.3 s (Extended Data Fig. 5c,d).

Mechanism in DI water

For the millimetre-scale mechanism in DI water (Extended Data Fig. 5a), under the applied voltage, electrolysis produces H+ near the anode and OH− near the cathode. Because the AAc-co-AAm network is pH-sensitive, acidic conditions convert –COO− to –COOH (reducing repulsion force and causing contraction), whereas alkaline conditions convert –COOH to –COO− (increasing repulsion and causing swelling). The bath-induced pH gradient therefore shrinks the anode-side region and swells the cathode-side region, resulting in bending towards the anode42.

By contrast, for the micrometre-scale mechanism in DI water (Extended Data Fig. 5b), the micrometre-scale hydrogel bends towards the cathode. Here fixed –COOH groups partially dissociate into –COO− and mobile H+. To maintain electroneutrality, the dissociated H+ ions largely remain confined within the hydrogel. Under an external field, these H+ ions migrate and accumulate on the cathode-facing side, locally lowering pH and inducing network contraction, which bends the hydrogel towards the cathode.

Mechanism in 0.15380 mol l−1 NaCl

For the millimetre-scale mechanism in 0.15380 mol l−1 NaCl (Extended Data Fig. 5c), in a saline environment, the external field drives the migration of all mobile ions in the bath. The fixed negative charges of the hydrogel influence ion partitioning at the gel–solution interface42,43,44,45, producing non-uniform ion concentrations across four regions (gel anode side: region 2; gel cathode side: region 1; solution anode side: region 4; solution cathode side: region 3). According to Flory’s theory51, the local osmotic pressure difference is \(\Delta {\rm{\pi }}={RT}\sum _{i}({c}_{i{\rm{g}}}-{c}_{i{\rm{s}}})\), in which Δπ is the pressure difference, cig is the ion concentration in the hydrogel, cis is the ion concentration in solution, R is the gas constant and T is temperature. Under steady-state conditions, the osmotic pressure difference on the anode side \(\Delta {{\rm{\pi }}}_{\text{Anode side}}={RT}\sum _{i}\,({c}_{i{\rm{region}}2}-{c}_{i{\rm{region}}4})\) exceeds that on the cathode side \(\Delta {{\rm{\pi }}}_{\text{Cathode side}}={RT}\sum _{i}({c}_{i{\rm{region}}1}-{c}_{i{\rm{region}}3})\). This imbalance causes the anode side to swell more, bending the hydrogel towards the cathode.

At the microscale, the dominant factor is the ion concentration gradient inside the hydrogel rather than at the interface. The relevant osmotic pressure is \(\Delta {{\rm{\pi }}}_{\text{Inside the hydrogel}}=RT\sum _{i}({c}_{i{\rm{r}}{\rm{e}}{\rm{g}}{\rm{i}}{\rm{o}}{\rm{n}}1}-{c}_{i{\rm{r}}{\rm{e}}{\rm{g}}{\rm{i}}{\rm{o}}{\rm{n}}2})\), which remains positive, making the cathode side swell more and driving bending towards the anode (Extended Data Fig. 5d).

Osmotic pressure analysis

Previously reported centimetre-scale and millimetre-scale hydrogel actuators operated in enclosed solution environments, in which the hydrogel was immersed in a bath equipped with two electrodes on the sidewalls44,50,52 (Extended Data Fig. 6a). In such systems, most ions in the bath are consumed in establishing concentration gradients across the hydrogel–solution interfaces, generating osmotic pressure differences that drive the hydrogel to bend towards the cathode42,43,44.

By contrast, our microscale hydrogel cilia operate within a localized region (200 µm × 200 µm × 90 µm) inside a much larger bath (4 cm × 4 cm × 3 mm) (Extended Data Fig. 6b). The large bath volume serves as an ion reservoir, allowing rapid diffusion of external ions into the actuation region (Extended Data Fig. 6b (ii)) and thereby preventing the formation of substantial ion concentration gradients or osmotic pressure differences across the hydrogel–solution interfaces.

To further validate this reason, osmotic-pressure simulations are performed under two distinct scenarios. Here osmotic pressure refers to the interfacial osmotic-pressure difference that governs centimetre-scale and millimetre-scale hydrogel actuation (as defined in Extended Data Fig. 5c). Direct comparison between simulations at different physical scales is not meaningful because the results would not be dimensionally consistent. Therefore, both simulations are conducted at the micrometre scale, with the only variable being the presence or absence of ion sources at the boundaries.

Case 1: macroscale actuation mimicked at the microscale (Extended Data Fig. 6c (i)). The simulation domain is a 2D region of 200 µm × 150 µm containing a hydrogel of dimensions 10 µm × 90 µm, reproducing the relative size ratio between the hydrogel and the bath in the centimetre–millimetre-scale experiments. No external ion source is applied, corresponding to the enclosed environment characteristic of macroscale systems.

Case 2: microscale actuation with ion exchange (Extended Data Fig. 6c (ii)). The geometry and parameters are identical to case 1, except that ion-source boundary conditions are imposed on three sides of the simulation domain, allowing continuous ion exchange between the system and the surroundings.

As shown in Extended Data Fig. 6d, introducing ion sources at the boundaries markedly reduced the osmotic-pressure difference across the hydrogel–solution interface. This result confirms that, compared with the macroscale condition, the osmotic contribution is greatly diminished at the microscale and is no longer the dominant mechanism governing hydrogel actuation.

Coupled electro-chemo-mechanical model

A fully coupled electro-chemo-mechanical model is developed here. The model simultaneously resolves: (1) ionic concentration distribution under an external electric field and the fixed charges of the hydrogel network; (2) the resulting osmotic body forces generated by non-uniform ion distributions; and (3) the gel deformation driven by these body forces. At this stage, the model focuses on explaining the deformation of hydrogels in ionic solutions under external electric fields and does not include fluid–structure interactions between the gel and the surrounding liquid.

Ionic concentration distribution

The Nernst–Planck equation describes the ion concentration in a charged hydrogel network environment:

$${J}_{\alpha }=-\,{D}_{\alpha }{\rm{\nabla }}{c}_{\alpha }-{z}_{\alpha }{\mu }_{\alpha }{c}_{\alpha }\nabla \O +{c}_{\alpha }v$$

(1)

in which Jα is the ion flux (mol (m2 s)−1), Dα is the diffusion coefficient of ion species α, cα is the concentration of ion species α, zα is the charge number of the ion (valence), μα is the mobility of the ion, Ø is the electric potential and v is the velocity of the fluid.

The change in ion concentration cα over time is governed by the continuity equation, which states that the change rate of the ion concentration is equal to the net flux of ions plus any sources or sinks (chemical reactions, for instance). It is expressed as:

$$\frac{{\rm{\partial }}{c}_{\alpha }}{{\rm{\partial }}t}+{\rm{\nabla }}\cdot {J}_{\alpha }={r}_{\alpha }({c}_{\beta })$$

(2)

in which \(\frac{{\rm{\partial }}{c}_{\alpha }}{{\rm{\partial }}t}\) is the ion concentration change rate, ∇·Jα represents the net flow of ions into or out of a region and rα(cβ) is a source term representing the creation or consumption of ion α owing to chemical reactions or other processes involving species β.

By substituting the Nernst–Planck equation (equation (1)) into the continuity equation (equation (2)), we obtain:

$$\frac{{\rm{\partial }}{c}_{\alpha }}{{\rm{\partial }}t}={\rm{\nabla }}\cdot [{D}_{\alpha }{\rm{\nabla }}{c}_{\alpha }+{z}_{\alpha }{\mu }_{\alpha }{c}_{\alpha }{\rm{\nabla }}\O -{c}_{\alpha }v]+{r}_{\alpha }({c}_{\beta })$$

(3)

The Poisson equation relates the electric potential distribution to the charge density in the system. It is given by

$${{\rm{\nabla }}}^{2}\O =-\frac{\rho }{\in }$$

(4)

in which ρ is the charge density (total charge per unit volume) and ε is the permittivity of the medium, expressed as εrε0, in which ε0 is the vacuum permittivity and εr is the relative permittivity of the medium.

In this work, the AAc-co-AAm hydrogel network contains –COOH groups, whose ionization equilibrium is influenced by H+ migration. This equilibrium affects the network charge and mobility of other ions. This equilibrium is included in the modelling by

$$\frac{{C}_{{\text{R-COO}}^{-}}\cdot {C}_{{{\rm{H}}}^{+}}}{{C}_{\text{R-COOH}}-{C}_{{\text{R-COO}}^{-}}}={K}_{a}=5.6\times {10}^{-5}$$

(5)

in which \({C}_{{\text{R-COO}}^{-}}\) is the concentration of the dissociated function group in the hydrogel, \({C}_{{{\rm{H}}}^{+}}\) is the H+ ion concentration in the hydrogel and CR-COOH is the concentration of the carboxyl group in the hydrogel. The value of CR-COOH is obtained from the initial hydrogel solution. Equations (3)–(5) collectively govern the ion concentration distribution.

Forces generated by non-uniform ion distributions

The force induced by ionic distributions can be separated into two contributions. The first contribution comes from the H+ effect (pH-dependent). H+ modulate the dissociation equilibrium of –COO− groups within the hydrogel network. Variations in the fixed –COO− concentration govern network swelling or contraction. The corresponding force can be expressed as:

$${f}_{{\rm{pH}}}=E(\,-\,{\rm{\nabla }}{({FC}}_{{\text{R-COO}}^{-}}))$$

(6)

in which F is the Faraday constant and E is the local electric field, which can be obtained from

$$E=-\,{\rm{\nabla }}\O .$$

(7)

The second contribution comes from the local pressure difference induced by ion concentration. According to Flory’s theory51, the local pressure relates to the ion concentration:

$${\rm{\pi }}={RT}\sum _{i}{c}_{i}$$

(8)

in which π is the local pressure, ci is the ion concentration, R is the gas constant and T is temperature.

The force generated by the non-uniform ion concentration can be expressed by the negative gradient of the local pressure as

$$f=-\,{\rm{\nabla }}{\rm{\pi }}=-\,{\rm{\nabla }}({RT}\sum _{i}{c}_{i}).$$

(9)

The local electric potential Ø, function group concentration \({C}_{{\text{R-COO}}^{-}}\) and ion concentration ci can be calculated from equations (3)–(5).

Mechanical deformation

The deformation of the hydrogel was modelled as a nonlinear hyperelastic material considering geometric nonlinearity. To describe its constitutive behaviour, the first-order compressible Ogden model was used. The strain energy density function is expressed as

$$W=\frac{\mu }{\alpha }({\lambda }_{1}^{\alpha }+{\lambda }_{2}^{\alpha }+{\lambda }_{3}^{\alpha }-3)+\frac{1}{D}{(J-1)}^{2}$$

(10)

in which λ1, λ2 and λ3 are the principal stretch ratios, J = λ1λ2λ3 is the volume ratio, μ and α are material constants and D is a compressibility parameter related to the bulk modulus K = 2/D. The Poisson’s ratio ν = 0.42 is chosen to define the degree of compressibility. The material parameters (μ, α) are determined by fitting the Ogden model to the data obtained from atomic force microscopy tests and are listed in Extended Data Table 2.

PIV analysis

The 2D in-plane velocity is calculated using the open-source software PIVlab53. Background subtraction is applied to the particle images to enhance image quality. The multipass fast Fourier transform window deformation algorithm is used to improve the accuracy of displacement estimation. The interrogation area is initially set to 64 × 64 pixels with a step size of 32 pixels, corresponding to a 50% overlap between adjacent interrogation windows. In the second pass, the interrogation area is reduced to 32 × 32 pixels with a step size of 16 pixels. For higher precision, sub-pixel displacements are estimated using a Gaussian 2 × 3 point estimator. The PIV flow pattern results are shown in Fig. 4.

It should be noted that the patterned electrodes on the substrate influence the accuracy of PIV analysis. For cases in Fig. 4a–c, in which the electrode density is relatively low, the PIV results are reliable. For case in Fig. 4d–g, the electrode density is higher and the trajectories of the PIV tracer particles overlap substantially with the underlying electrodes, leading to larger errors in velocity quantification.

2D flow simulations

For the numerical simulation of the interaction between the cilia array and fluid, we use the hybrid finite difference/finite element immersed boundary method54, implemented in the open-source software IBAMR, a widely tested C++ framework for the immersed boundary method. The immersed boundary formulation of the problem describes the momentum and velocity of the coupled fluid–structure system in Eulerian form, whereas the deformation and elastic response of the immersed structure is in the Lagrangian form. This study uses 2D simulations, as they effectively capture the flow patterns. Let \({\bf{x}}=({x}_{1},{x}_{2})\in \Omega \subset {{\mathbb{R}}}^{2}\) denote Cartesian physical coordinates, in which Ω represents the physical region that is occupied by the coupled fluid–structure system, let \({\bf{X}}=({X}_{1},{X}_{2})\in W\subset {{\mathbb{R}}}^{2}\) denote Lagrangian material coordinates that are attached to the structure, in which W is the Lagrangian domain, and let χ(X, t) ∈ Ω denote the physical position of material point X at time t. The strong form of the equations of motion is:

$$\begin{array}{l}\rho \frac{{\rm{D}}u}{{\rm{D}}t}({\bf{x}},t)=-\,{\rm{\nabla }}p({\bf{x}},t)+\mu {{\rm{\nabla }}}^{2}{\bf{u}}({\bf{x}},t)+{{\bf{f}}}^{{\bf{c}}}({\bf{x}},t)\\ {\rm{\nabla }}\cdot {\bf{u}}({\bf{x}},t)=0\end{array}$$

(11)

$$\begin{array}{l}\begin{array}{c}{{\bf{f}}}^{{\bf{c}}}({\bf{x}},t)={\int }_{U}{{\rm{\nabla }}}_{{\bf{X}}}\cdot {{\mathbb{P}}}^{{\bf{e}}}({\bf{X}},t)\delta ({\bf{x}}-\chi ({\bf{X}},t)){{\bf{f}}}^{{\bf{e}}}({\bf{x}},t){\rm{d}}{\bf{X}}\\ \,-{\int }_{{\rm{\partial }}U}{{\mathbb{P}}}^{{\bf{e}}}({\bf{X}},t){\bf{N}}({\bf{X}})\delta ({\bf{x}}-\chi ({\bf{X}},t)){{\bf{f}}}^{{\bf{e}}}({\bf{x}},t){\rm{d}}{\bf{X}}\\ \,\frac{{\rm{\partial }}\chi }{{\rm{\partial }}t}({\bf{X}},t)={\int }_{\Omega }{\bf{u}}({\bf{x}},t)\delta ({\bf{x}}-\chi ({\bf{X}},t)){\rm{d}}{\bf{X}}\end{array}\end{array}$$

(12)

in which ρ is the mass density, u(x, t) is the Eulerian velocity field, μ is the dynamic viscosity, fe(x, t) is the Eulerian elastic force density, \({{\mathbb{P}}}^{{\bf{e}}}({\bf{X}},t)\) is the first Piola–Kirchhoff elastic stress tensor, δ(x) is the 2D delta function and N(X) is the normal vector along the fluid–solid interface. In the computations of this study, the physical domain is Ω = [−L, L][−L, L], in which L is 600 μm for simulations corresponding to fluid experiments 1–3 and 400 µm for fluid experiments 4–7. A zero-gradient boundary condition is applied to boundaries. A staggered-grid finite difference scheme is used to discretize the incompressible Navier–Stokes equations in space. The spatial resolution is Δx = L/128 for simulations 1–3 and Δx = L/64 for simulations 4–7. The total number of the Cartesian grid is \({\mathcal{O}}({10}^{5})\). The circular section of the cilia is discretized into a mesh of triangular elements with an average node spacing of ΔX = L/128. Time-stepping is performed using an implicit scheme proposed by Newren et al.55. The time step size is adjusted to satisfy the Courant–Friedrichs–Lewy condition, with a stability number of approximately 0.1.

In these 2D flow-field simulations, solid spheres with a diameter of 10 µm are used to approximate the hydrogel cilia. Under the prescribed motions, interactions between the sphere edges and the surrounding fluid generated relatively high flow velocities. However, in reality, hydrogel cilia are porous and soft materials and their interaction with the fluid under identical motions does not produce flow velocities of the same magnitude as in the simulations. The primary aim of the simulation is to qualitatively predict flow-field patterns, not to precisely quantify flow velocities, so the simulated flow patterns are shown in this work (Fig. 4a–g (2)).