1. Introduction
Ciliate-like structures are abundant in nature, manifesting in various forms such as the cilia found on plant leaves, paramecia, coral surface structures (Eloy & Lauga Reference Eloy and Lauga2012; Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky Winter, Vardi and Stocker2014; Bouhouche et al. Reference Bouhouche, Valentine, Le Borgne, Lemullois, Yano, Lodh, Nabi, Tassin and Van Houten2022), as well as in the respiratory system and oviducts (Koyama, Shi & Fujimori Reference Koyama, Shi and Fujimori2019; Causa et al. Reference Causa, Das, Feriani, Kotar and Cicuta2025). These structures interact with the surrounding fluid, impacting the morphology of the flow field and playing crucial roles in diverse biological activities (Brennen & Winet Reference Brennen and Winet1977; Kim et al. Reference Kim, Bird, Van Parys, Breuer and Powers2003; Liu, Powers & Breuer Reference Liu, Powers and Breuer2011, Reference Liu, Costello and Kanso2025a , Reference Liu, Man, Costello and Kansob ). For instance, the active oscillation of sperm cilia effectively aids in their own movement (Gaffney et al. Reference Gaffney, Gadêlha, Smith, Blake and Kirkman Brown2011; Lauga & Goldstein Reference Lauga and Goldstein2012). Furthermore, the oscillating cilia found on paramecia and coral surfaces generate vortices, enabling the exchange of nutrients and gases with the environment in a captivating display of fluid dynamics. Behind these phenomena, the collective role of cilia plays an important role in the life activities of these organisms. Cilia, found in diverse biological structures like the respiratory system and oviducts, employ coordinated oscillations to achieve tasks as varied as dust and egg transportation, fluid mixing and nutrient exchange, often characterised as metachronal waves (Gilpin, Bull & Prakash Reference Gilpin, Bull and Prakash2020; Mesdjian et al. Reference Mesdjian, Wang, Gsell, D’ortona, Favier, Viallat and Loiseau2022; von Kenne et al. Reference von Kenne, Schmelter, Stark and Bär2025).
By harnessing the power of metachronal waves in artificial cilium systems (Coq et al. Reference Coq, Bricard, Delapierre, Malaquin, Du Roure, Fermigier and Bartolo2011; Elgeti & Gompper Reference Elgeti and Gompper2013; Dong et al. Reference Dong, Lum, Hu, Zhang, Ren, Onck and Sitti2020; Gu et al. Reference Gu2020), one can observe remarkable potential to enhance material transport and fluid dynamics. Similarly, paramecia and ciliated structures on coral surfaces transport oxygen and nutrients by generating vortices (Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky Winter, Vardi and Stocker2014; Bouhouche et al. Reference Bouhouche, Valentine, Le Borgne, Lemullois, Yano, Lodh, Nabi, Tassin and Van Houten2022). Studies reveal that longitudinal waves maximise mucus transport, while transverse waves minimise energy dissipation by cilia (Mesdjian et al. Reference Mesdjian, Wang, Gsell, D’ortona, Favier, Viallat and Loiseau2022). An important theoretical framework for ciliary dynamics is the scallop theorem (Purcell Reference Purcell1977; Lauga Reference Lauga2011b ), which asserts that reciprocal motion of rigid bodies in low Reynolds number flows, such as cilia, cannot generate net fluid transport due to the time-reversible nature of the flow (Qiu et al. Reference Qiu, Lee, Mark, Morozov, Münster, Mierka, Turek, Leshansky and Fischer2014; Esparza López & Lauga Reference Esparza López and Lauga2023). As a result, ciliary systems rely on asymmetric flapping motions, like metachronal waves, to efficiently drive or pump fluid (Hall & Clarke Reference Hall and Clarke2020; Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2024; Xia et al. Reference Xia, Hu, Wei, Chen, Peng and Yang2024). Nevertheless, many cilia operate within chaotic background flows, particularly in artificial cilium systems at relatively larger scales. These chaotic flows introduce additional complexities that challenge traditional models of the ciliary dynamics. For example, enhanced diffusion can be seen for a swimmer in a fluctuating environment (Lauga Reference Lauga2011a ). However, the precise coupling between cilium carpet and chaotic flows remains poorly understood, and how cilia strategically harness oscillation energy to modify flow structures for enhanced transport efficiency is still unknown, where there is a modified interaction between the flexible structure and the flow (Médéric & Mahadevan Reference Médéric and Mahadevan2005; Shelley & Zhang Reference Shelley and Zhang2011; Alben Reference Alben2015; Guttag & Reis Reference Guttag and Reis2017; Ganedi et al. Reference Ganedi, Oza, Shelley and Ristroph2018; Mathai et al. Reference Mathai, Das, Naylor and Breuer2023). Recent progress on flexible fibres in turbulence – reviewed in Marchioli, Rosti & Verhille (Reference Marchioli, Rosti and Verhille2025) – highlights how flexibility enables local coupling with multiscale flow structures. In this context, Mathai et al. (Reference Mathai, Tzezana, Das and Breuer2022) demonstrated that a compliant membrane hydrofoil passively cambered, stabilised leading-edge vortices, and achieved up to 160 % higher power than rigid foils. Building on this, Upfal et al. (Reference Upfal, Zhu, Handy-Cardenas and Breuer2025) identified constant- and dynamic-camber regimes in which elasticity and extensibility stabilise vortices and soften stall. These insights provide a concise framework for interpreting elastic filament–flow interactions in chaotic convection.
Strikingly, here, we address the open question of how cilium-synchronised oscillation can be utilised in the chaotic flow to optimise scalar transport using thermal systems as an example. Traditionally, the laminar boundary layer imposes a limitation on scalar transport; however, we demonstrate that oscillating cilia can induce resonance within the boundary layer, which enhances transport at intermediate frequencies. Our results demonstrate a mechanistic transition: at low to medium frequencies, cilia uplift fluid, augmenting heat transfer through vertical displacement. Beyond a critical frequency, an intensified flow dynamics generates vortices from unstable shear layers, drastically improving fluid mixing and thermal efficiency. Our findings demonstrate how cilia can exploit chaotic background flows to enhance scalar transport.
To systematically investigate these effects, we consider a paradigm of buoyancy-driven flow – Rayleigh–Bénard (RB) convection, i.e. convective flow inside a cubic box heated from below and cooled from above (Siggia Reference Siggia1994; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Shishkina Reference Shishkina2021; Lohse & Shishkina Reference Lohse and Shishkina2023, Reference Lohse and Shishkina2024; Xia et al. Reference Xia, Chong, Ding and Zhang2025). Previous studies have explored various strategies to regulate heat transport in thermal convection (Wang et al. Reference Wang, Mathai and Sun2019, Reference Wang, Zhou and Sun2020; Jiang et al. Reference Jiang, Zhu, Wang, Huisman and Sun2020, Reference Jiang, Wang, Liu and Sun2022; Xia et al. Reference Xia, Huang, Xie and Zhang2023; Xu et al. Reference Xu, Zhang, Brandt, Song and Shishkina2025). Motivated by this, we introduce cilia uniformly distributed on the hot plate, with spacing
$\varLambda$
and length
$\mathcal{L}$
comparable to the thermal boundary-layer thickness (see figure 1
$a$
). In this work, we use ‘cilia’ as a shorthand for wall-mounted oscillating filament arrays, encompassing both biological cilia and engineered flexible fibres. The paper is organised as follows: The numerical set-ups are described in § 2, followed by the heat transfer properties, flow structure and underlying mechanism in § 3. The conclusions are presented in § 4.
2. Numerical method
In the present study, we consider the flow within a two-dimensional square box and perform direct numerical simulations (DNS) of thermal convection with an oscillating filament array fixed at the bottom plate. The non-dimensional governing equations under the Boussinesq approximation are
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \partial _t \boldsymbol {u} + (\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u} = -\boldsymbol{\nabla }p + \sqrt {\frac {\textit{Pr}}{\textit{Ra}}} {\nabla} ^2 \boldsymbol {u} + \theta \boldsymbol {e}_z + \boldsymbol {f}, \\[-28pt] \nonumber \end{align}
\begin{align} \partial _t \theta + (\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \theta = \frac {1}{\sqrt {\textit{RaPr}}} {\nabla} ^2 \theta , \\[-2pt] \nonumber \end{align}
where
$t$
is the dimensionless time scaled by the free-fall time
$\sqrt {H / (\alpha g \Delta )}$
,
$\boldsymbol {u}$
is the fluid velocity scaled by the free-fall velocity
$\sqrt {\alpha g \Delta H}$
and
$\theta$
is the temperature scaled by the imposed temperature difference
$\Delta$
;
$\boldsymbol {e}_z$
is the unit vector along the vertical direction, and
$\boldsymbol {f}$
is the feedback force exerted by the filament on the fluid. All length scales are normalised by the cell height
$H$
. The Rayleigh number
$\textit{Ra}$
and Prandtl number
$\textit{Pr}$
are defined as
\begin{align} \textit{Ra} = \frac {\alpha g \Delta H^3}{\nu \kappa }, \qquad \textit{Pr} = \frac {\nu }{\kappa }, \end{align}
where
$\alpha$
is the thermal expansion coefficient,
$g$
the gravitational acceleration,
$\nu$
the kinematic viscosity and
$\kappa$
the thermal diffusivity of the fluid. At the boundaries, no-slip velocity conditions are applied at all walls, with constant temperatures
$\theta _{\textit{top}} = 0$
and
$\theta _{\textit{bot}} = 1$
at the top and bottom plates, respectively. Adiabatic conditions are applied at the sidewalls.
For the filament, we model inextensible filaments represented by Lagrangian nodes. The dynamics of a thin, flexible filament is governed by (2.5) and (2.6) under the inextensibility constraint (Huang, Shin & Sung Reference Huang, Shin and Sung2007; Pinelli et al. Reference Pinelli, Omidyeganeh, Brucker, Revell, Sarkar and Alinovi2017; Banaei, Rosti & Brandt Reference Banaei, Rosti and Brandt2020). The dimensionless equations and the inextensibility condition are
\begin{align} \frac {\partial ^{2} \boldsymbol{X}}{\partial t^{2}}=\frac {\partial ^{2} \boldsymbol{X_{f}}}{\partial t^{2}}+\frac {\partial }{\partial s} \left(T\frac {\partial \boldsymbol{X}}{\partial s} \right)-\frac {\partial ^{2}}{\partial s^{2}} \left(B\frac {\partial ^{2} \boldsymbol{X}}{\partial s^{2}} \right)-\boldsymbol{F_{\kern-1.5pt f}}+\boldsymbol{F_{c}}, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol{X}}{\partial s} \boldsymbol{\cdot }\frac {\partial \boldsymbol{X}}{\partial s} = 1, \\[-2pt] \nonumber \end{align}
where
$\boldsymbol{X}$
is the filament position,
$s$
denotes arclength,
$ {\partial \boldsymbol{X_{\kern-1.5pt f}}}/{\partial t}$
is the fluid velocity at the filament node,
$T$
is the dimensionless tension force,
$B$
is the dimensionless bending rigidity,
$\boldsymbol{F_{\kern-1.5pt f}}$
is the fluid force on the filament and
$\boldsymbol{F_c}$
is the repulsive force between filament nodes to prevent overlap and self-entanglement. These equations, along with the inextensibility condition, are also described in detail by Huang et al. (Reference Huang, Shin and Sung2007) and Banaei et al. (Reference Banaei, Rosti and Brandt2020).
We induce oscillations in the cilia by changing the boundary conditions and investigate how this modulation affects global heat transport in the system. The boundary condition at the root of cilia is given by
\begin{align} \frac {\partial \boldsymbol{X}}{\partial s} = \left[\sin \left (\frac {\pi }{2} \sin (2 \pi \varOmega t)\right ), \cos \left (\frac {\pi }{2} \sin (2 \pi \varOmega t)\right ) \right]\!, \end{align}
where
$\varOmega$
is the oscillation frequency. For thermal coupling, we assume rapid heat exchange between each filament and the surrounding fluid; accordingly, the filament temperature equals the local fluid temperature (the filament is locally isothermal).
For the cilium motion and arrangement, we introduce five additional dimensionless control parameters: the dimensionless tension
$T = \mathcal{T} / (\rho _f \alpha g \Delta H)$
, the dimensionless bending rigidity
$B = \mathcal{B} / (\rho _f \alpha g \Delta H^3)$
, the dimensionless frequency
$\omega = \sqrt {\varOmega ^2 H / (\alpha g \Delta )}$
, the dimensionless length
$L = \mathcal{L} / H$
and the dimensionless spacing
$\lambda = \varLambda / H$
, where
$\rho _f$
is the filament linear density (mass per unit length) (Banaei et al. Reference Banaei, Rosti and Brandt2020),
$\mathcal{T}$
is the cilium tension and
$\mathcal{B}$
is the cilium bending rigidity. Unless stated otherwise, we assume neutral buoyancy (the filament density matches that of the surrounding bulk fluid) and fix the bending-rigidity parameter at
$B=0.2$
.
(
$a$
) Schematic illustrating enhanced heat transport in thermal convection through active cilia, promoting efficient thermal energy redistribution between the heating and cooling plates. (
$b$
) Normalised Nusselt number
$ \textit{Nu}(\omega )/ \textit{Nu}_{\textit{RB}}$
as a function of oscillation frequency
$\omega$
for various cilium lengths
$L$
at
$\textit{Ra}$
=
$10^7$
and
$\textit{Pr}$
= 4.38. Here,
$ \textit{Nu}_{\textit{RB}}$
denotes the Nusselt number of standard RB convection (
$ \textit{Nu}_{\textit{RB}}$
= 13.05 for two dimensions, 16.8 for three dimensions). (
$c$
–
$e$
) Behaviour of plumes within the cilium layer at various oscillation frequencies:
$\omega = 0.05$
in (
$c$
),
$\omega = 1$
in (
$d$
),
$\omega = 40$
in (
$e$
). The middle panels show the vertical slice of temperature at the mid-length of the cilia for each frequency, while the upper panels depict the vertical slice of
$w\theta$
at the same location.
To solve the governing equations, we use the finite-difference solver AFiD (Verzicco & Orlandi Reference Verzicco and Orlandi1996; Ostilla-Monico et al. Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015) with a high-performance message passing interface. The equations are solved using a multi-scalar second-order finite-difference method coupled with a fractional third-order Runge–Kutta scheme. Fluid–filament coupling is achieved via the immersed boundary method. The volume force exerted by the filaments on the fluid is computed by convolving the singular forces from the Lagrangian nodes onto the Eulerian mesh (Huang et al. Reference Huang, Shin and Sung2007; Banaei et al. Reference Banaei, Rosti and Brandt2020). Our code has been used to simulate various single-phase and multiphase turbulent flow problems (Guo et al. Reference Guo, Wu, Wang, Zhou and Chong2023; Zhao et al. Reference Zhao, Zhao, Zhou and Chong2024). For fibre-dynamics validation, our results agree with Huang et al. (Reference Huang, Shin and Sung2007) (details in Appendix A). We performed DNS of RB convection with cilium oscillations in both two-dimensional and three-dimensional square boxes. We fixed
$\textit{Ra}=1 \times 10^{7}$
,
$\textit{Pr}=4.38$
(corresponding to water conditions) and
$\lambda =0.1$
. A series of DNS runs is conducted over the parameter ranges
$0.05 \leqslant \omega \leqslant 80$
and
$0.04 \leqslant L \leqslant 0.06$
. For the three-dimensional cases, we set the aspect ratio of the box to 0.15 for width to height, with a length-to-height ratio of 1 : 1. Cilia are arranged with 30 fibres along the x-direction, while the distribution in the y-direction follows the same spacing as in the two-dimensional cases.
3. Results
Figure 1(
$b$
) depicts the variation of the global heat flux
$ \textit{Nu}$
with the cilium oscillation frequency
$\omega$
at
$\textit{Ra}=10^7$
and
$\textit{Pr} = 4.38$
. It is shown that the dependence of
$ \textit{Nu}$
on frequency follows a similar trend, with the
$\omega$
-dependence of
$ \textit{Nu}$
showing three distinct regimes. In the low-frequency range, the
$ \textit{Nu}$
value is slightly lower than that of classical RB convection, with oscillating cilia having minimal impact on heat transport. As the frequency
$\omega$
increases, the
$ \textit{Nu}$
value begins to increase, indicating a substantial enhancement in heat transport due to the oscillating cilia. In particular, in the medium-frequency regime, the
$ \textit{Nu}$
-trend is not monotonic as
$\omega$
increases. A local maximal heat flux is observed, indicating the emergence of a local optimal frequency
$\omega _{{c}}$
at which the peak value
$ \textit{Nu}_{{c}}$
is located. Notably, compared with the RB convection, a remarkable
$ \textit{Nu}$
-enhancement rate reaches up to
$50\,\%$
. Moreover, we compare the trend of
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
for different lengths (
$L$
) of the cilia. One can see that the onset frequency is sensitive to
$L$
, and the regime transition shifts toward lower frequencies as the length of the cilia increases.
To gain more insight into the effect of oscillating cilia, we analyse how the flow structures within the ciliary layer are modified by the cilium motions. Figure 1(
$c$
–
$e$
) presents the instantaneous temperature fields within the ciliary layer for frequencies of
$\omega =0.05$
in the low-frequency regime,
$\omega =1$
in the medium-frequency regime and
$\omega =40$
in the high-frequency regime, all at a fixed spacing of
$\lambda =0.1$
. In the low-frequency regime (see figure 1
$c$
), the oscillating cilia have minimal impact on plume emissions (see Supplementary Materials Movie S1 are available at https://doi.org/10.1017/jfm.2025.11000). However, when the frequency exceeds the critical frequency in the medium-frequency regime (see figure 1
$d$
), the oscillation of the cilia disrupts the thermal boundary layer and triggers substantial plume emissions (see Supplementary Materials Movie S2 for details). Under the shear effect of the cilia, the plumes oscillate along with the cilia; as the frequency further increases in the high-frequency regime, the plumes cease to oscillate and instead stabilise with reduced amplitude (see figure 1
$e$
and Supplementary Materials Movie S3 ). What physics then governs the underlying mechanism of cilium modulation of the plume dynamics?
(
$a$
–
$e$
) Morphological evolution of cilia and plumes over one oscillation cycle: active cilia sweep through the temperature field, uplifting the thermal boundary layer and generating plumes. (
$f$
) Detailed representation of the strong shear effect resulting in a localised uplift mechanism in the thermal boundary layer during the cilia’s interactions with the temperature field. (
$g$
) Normalised Nusselt number
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
vs
$\omega /\omega _{{c}}$
for various cilium lengths
$L$
. (
$h$
)Normalised Reynolds number
$ \textit{Re}(\omega )/ \textit{Re}_{\textit{RB}}$
vs
$\omega /\omega _{{c}}$
for various cilium lengths
$L$
. The vertical dashed line represents
$\omega /\omega _{{c}}$
= 1 indicates the onset of enhancement in both
$ \textit{Nu}$
and
$ \textit{Re}$
. Here,
$ \textit{Re}=\sqrt {\langle u^2+w^2\rangle }H / \nu$
, with
$\langle \boldsymbol{\cdot }\rangle$
denoting spatial and temporal averaging, and
$ \textit{Re}_{\textit{RB}}$
is the Reynolds number for standard RB convection (
$ \textit{Re}_{\textit{RB}}$
= 184.7 for two dimensions, 85.7 for three-dimensions).
First, we focus on how oscillating cilia disrupt the thermal boundary layer, thereby triggering significant plume emissions. Figure 2(
$a$
–
$e$
) illustrates the morphological evolution of cilia and plumes over one oscillation cycle: active cilia sweep through the temperature field, uplifting the thermal boundary layer and generating plumes, which suggests a complex interaction between the cilium oscillation and the flow dynamics within the boundary layer (Lee et al. Reference Lee, Oh, Chung, Chang, Choi, Yagawa and Kim2011; Zhu, He & Zhang Reference Zhu, He and Zhang2014; Stein et al. Reference Stein, De Canio, Lauga, Shelley and Goldstein2021). As shown in figure 2(
$f$
), when cilia sweep through the thermal boundary layer, the shear forces exerted on the surrounding fluid can uplift and destabilise it. Specifically, at low oscillation frequencies, the sweeping time scale
$\tau _{{os}}\sim ( {\delta _{\textit{th}}}/{\omega L})$
exceeds the viscous dissipation time scale
$\tau _{\textit{dis}}\sim ( {\delta _{\textit{th}}^2}/{\nu })$
, resulting in rapid dissipation of fluid disturbances without significant disruptions to the flow or temperature fields. In contrast, as the frequency increases and the sweeping time approaches the dissipation time scale, the shear-induced disturbances are not fully dissipated, leading to the continuous acceleration of fluid microstructures. Thus, by balancing these two time scales, we can derive the critical condition
$\omega _{{c}}\sim ( {\nu }/{\delta _{\textit{th}}L})$
. Here, the thermal boundary-layer thickness is obtained from the global Nusselt number as
$\delta _{\textit{th}}=1/(2\textit{Nu})$
.
Next, to confirm this proposed physical picture of the onset of
$ \textit{Nu}$
enhancement, we replot
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
as a function of the normalised frequency
$\omega /\omega _{{c}}$
for various values of
$L$
in figure 2(
$g$
). The figure shows that the normalised
$ \textit{Nu}$
does indeed begin to increase from 1 when
$\omega \gt \omega _{{c}}$
. This trend is robust across a broad parameter range: additional simulations spanning Rayleigh number
$\textit{Ra}$
, filament rigidity
$B$
and inter-cilium spacing
$\lambda$
exhibit the same qualitative dependence of
$ \textit{Nu}$
on
$\omega /\omega _{{c}}$
(see Appendix B). In addition, we present a plot of the normalised Reynolds number
$ \textit{Re}(\omega )/\textit{Re}_{\textit{RB}}$
versus the oscillation frequency
$\omega$
for various cilium lengths
$L$
in figure 2(
$h$
). The qualitative trends indicate that
$ \textit{Re}$
exhibits minimal dependence on frequency at low oscillation frequencies, while the Reynolds number remains close to that observed in standard RB convection. However, when the frequency exceeds a critical value,
$ \textit{Re}$
increases significantly. Furthermore, the enhancement of
$ \textit{Re}$
can also be observed for conditions where
$\omega /\omega _{{c}}\gt$
1, reinforcing that our explanation for the onset behaviour is applicable in the
$ \textit{Re}$
regime as well (see Appendix B for additional parameters).
It is noteworthy that, during the initial stage of enhancement, the value of
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
exhibits a similar trend for different values of
$L$
. However, with further increases in the ratio
$\omega /\omega _{{c}}$
, a second transition point is observed (see figure 2
$g$
). For example, at
$L=0.6$
in figure 2(
$g$
), the normalised
$ \textit{Nu}$
begins to rise above 1 when
$\omega \gt \omega _{{c}}$
. Nevertheless, as the ratio
$\omega /\omega _{{c}}$
continues to increase (for
$\omega /\omega _{{c}}\gt 10$
), the normalised
$ \textit{Nu}$
starts to level off and even exhibits a slight decrease before increasing again. This behaviour is also reflected in the
$ \textit{Re}$
versus
$\omega$
relationship shown in figure 2(
$h$
).
In order to elucidate the transition from the medium-frequency regime to the high-frequency regime (the second transition), we first investigate the influence of ciliary oscillations on heat flux transport within the boundary layer, represented by the parameter
$\sqrt {\textit{RaPr}}\langle w\theta \rangle _{L}$
. As shown in figure 3(
$a$
), the time series of spatially averaged heat flux
$\sqrt {\textit{RaPr}}\langle w\theta \rangle _{L}$
within the ciliary layer over five cycles reveals a distinct behaviour: in the medium-frequency regime, characterised by
$\omega =1$
, each cycle displays two maxima, while in the high-frequency regime (with
$\omega =40$
), only a single maximum is observed. This finding is further confirmed by their corresponding spectral analyses; at
$\omega =1$
, the dominant frequency is
$f=2\omega$
, while at
$\omega =40$
, it is
$f=\omega$
(see figure 3
$c$
,
$d$
). To gain deeper insights into the flow mechanisms associated with these two dominant frequencies, we present the temperature and velocity fields at the peak moments of thermal transport over one cycle, as shown in figure 3(
$e$
,
$f$
). At
$\omega =1$
, the two peaks correspond to the occurrence of two significant uplift events induced by ciliary oscillations (see figure 3
$e$
), as evidenced by the temperature profiles in the midheight of the cilia presented in figure 1(
$d$
). In contrast, at
$\omega =40$
, heat transport is predominantly driven by localised vortices induced by high-frequency oscillations, rather than uplift phenomena associated with ciliary motion. The oscillatory behaviour of the cilia facilitates the emergence of shear vortices during the transition between the first and second half-cycles. This effect can be attributed to the stable shear flow established by the cilia during the initial phase of the cycle; the subsequent rapid directional shift in ciliary motion generates opposing shear forces, thereby promoting the development of localised vortices (see figure 3
$f$
). Ultimately, this observation underscores the transition between uplift-driven and vortex-driven heat transfer mechanisms.
(
$a$
,
$b$
) Time series of
$\sqrt {\textit{RaPr}}\langle w\theta \rangle _{L}$
for different frequencies:
$\omega =1$
(
$a$
) and
$\omega =40$
(
$b$
). Here,
$\langle \boldsymbol{\cdot }\rangle _{L}$
denotes averaging over the region within the filament height. Panels (
$c$
) and (
$d$
) show the corresponding spectral plots for (
$a$
) and (
$b$
), respectively. The black dashed line indicates the proportions at the frequencies
$f=\omega$
and
$f=2\omega$
. Panels (
$e$
) and (
$f$
) present the temperature field (colour) and velocity field (arrows) within the ciliary layer during the cilium oscillation period at moderate frequency (
$\omega =1$
) and high frequency (
$\omega =40$
). Panel(
$e$
) depicts the event of lifting twice via ciliary oscillation at moderate frequency, while panel(
$f$
) illustrates the event of rolling vortices via ciliary oscillation during the cilium oscillation period at high frequency. (
$g$
) Normalised Reynolds number
$ \textit{Re}_{l}(\omega )/\textit{Re}_{{s}}$
vs
$\omega$
for various cilium lengths
$L$
. Here,
$ \textit{Re}_{l}=\sqrt {\langle u^2+w^2\rangle _{l}}H / \nu$
is the local Reynolds number, where
$\langle \boldsymbol{\cdot }\rangle _{l}$
denotes averaging over the cilium region and in time. (
$h$
,
$i$
) Normalised Reynolds number
$ \textit{Re}(\omega )/\textit{Re}_{\textit{RB}}$
and Nusselt number
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
vs
$\omega /\omega _{{c2}}$
for different cilium lengths
$L$
.
We further investigate whether the transition from uplift-driven to vortex-driven heat transfer mechanisms arises from the increasing dominance of shear flow dynamics within the ciliary layer. Figure 3(
$g$
) shows the Reynolds number within the ciliary layer, normalised by the shear Reynolds number
$ \textit{Re}_{{l}}/\textit{Re}_{{s}}$
, as a function of frequency
$\omega$
. The local Reynolds number is
$ \textit{Re}_{l} = \sqrt {\langle u^2 + w^2 \rangle _{l}} H / \nu$
, with
$\langle \boldsymbol{\cdot }\rangle _{l}$
denoting averaging over the ciliary region. The flow intensity driven by the cilia is characterised by the shear Reynolds number
$ \textit{Re}_{{s}} = \omega \textit{LH} / \nu$
. At low frequencies, flow within the ciliary layer is dominated by large-scale circulation (classical RB convection). As frequency increases, shear from ciliary oscillations becomes dominant. Resonance occurs when the shear Reynolds number approaches the local Reynolds number
$ \textit{Re}_{{l}}$
, at
$ \textit{Re}_{{l}} / \textit{Re}_{{s}} = 1$
. From this, we derive the critical frequency
$\omega _{{c2}}=\textit{Re}_{{l}}\nu /\textit{LH}$
, and we replot
$ \textit{Re}(\omega )/\textit{Re}_{\textit{RB}}$
and
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
as a function of the normalised frequency
$\omega /\omega _{{c2}}$
for various values of
$L$
in figure 3(
$h$
,
$i$
). The figure demonstrates that both the Reynolds number
$ \textit{Re}$
and normalised Nusselt number
$ \textit{Nu}$
begin to transition from uplift-driven to vortex-driven heat transfer when
$\omega \gt \omega _{{c2}}$
. This also effectively anticipates the notable reversal in the trend depicted by the final data point in figures 3(
$h$
) and 3(
$i$
).
4. Conclusion
In summary, we introduce a novel approach to enhance scalar transport by harnessing the oscillation energy of cilia in chaotic flows, demonstrated via DNS in both two and three dimensions. Our result demonstrates that periodic ciliary oscillations, tuned to specific frequencies, significantly boost scalar transport by overcoming the constraints of laminar boundary layers. The enhanced transport is primarily due to boundary-layer resonance and a transition in the heat transport mode. At medium to low frequencies, the cilia effectively uplift the fluid, improving heat transport through direct vertical displacement. As the frequency increases, the ability of the cilia to induce uplift diminishes, leading to localised vortex formation. This shift is driven by the increasing dominance of the shear flow dynamics within the ciliary layer. Once the oscillation frequency reaches a critical threshold, the flow intensity becomes strong enough to develop vortices, generated by the unstable shear flow from the rapid directional changes of the oscillating cilia. These vortices significantly improve fluid mixing and heat transport efficiency. Thus, it is the intensified flow within the ciliary layer, driven by high-frequency oscillations, that triggers the transition to a vortex-dominated transfer mechanism. Notably, we employ a symmetric sinusoidal stroke; non-reciprocal strokes should intensify power-phase sweeping, enhance transport and slightly shift the optimal frequency while preserving the boundary-layer sweeping/diffusive scaling. This concept can be similarly applied to other chaotic systems where another force agent drives the chaotic flow rather than buoyancy.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.11000.
Funding
This work is supported by NSFC Excellence Research Group Program for ‘Multiscale Problems in Nonlinear Mechanics’ (No.12588201), and it was also supported by the Natural Science Foundation of China under grant nos. 12372219, 52450304, 12432011, 12422208, 12421002, 12372220 and 12502255, the China Postdoctoral Science Foundation under Grant No. 2024M761631 and the New Cornerstone Science Foundation through the New Cornerstone Investigator Program and the XPLORER PRIZE.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation of fibre-dynamics numerical solution
We validate our solver by simulating the oscillation of a filament under gravity in a quiescent environment, as described in Huang et al. (Reference Huang, Shin and Sung2007). For finite bending stiffness (
$B = 0.01$
), the filament oscillates at its natural frequency. As shown in figure 4, our results agree well with Huang et al. (Reference Huang, Shin and Sung2007), using a resolution of 100 Lagrangian points per unit length.
(
$a$
) Superposition of filament positions at successive times over half an oscillation period. (
$b$
) Time history of the free end position of the filament, with solid dots taken from Huang et al. (Reference Huang, Shin and Sung2007).
(
$a$
–
$c$
) Normalised Nusselt number
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
vs oscillation frequency
$\omega /\omega _c$
for various
$\textit{Ra}$
(
$a$
), spacing
$\lambda$
(
$b$
) and filament rigidity
$B$
(
$c$
). Panels (
$d$
–
$f$
) show
$ \textit{Re}/\textit{Re}_{\textit{RB}}$
vs
$\omega /\omega _c$
for different
$\textit{Ra}$
(
$d$
), spacing
$\lambda$
(
$e$
) and filament rigidity
$B$
(
$f$
). Here,
$ \textit{Nu}_{\textit{RB}}$
is the Nusselt number for standard RB convection. All data are for two-dimensional simulations.
Appendix B. Critical transition for different parameters
Figure 5 shows
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
and
$ \textit{Re}(\omega )/\textit{Re}_{\textit{RB}}$
versus the normalised frequency
$\omega /\omega _{{c}}$
for various parameters. At low
$\omega$
,
$ \textit{Nu}(\omega )/\textit{Nu}_{\textit{RB}}$
and
$ \textit{Re}(\omega )/\textit{Re}_{\textit{RB}}$
remain approximately constant, and then increase sharply beyond a critical frequency (
$\omega /\omega _{{c}} \gt 1$
), supporting the general applicability of our explanation. It is worth noting that the current results were obtained within a limited range, and extending this range could affect the heat transport magnitude. The simulations focus on relatively stiff filaments, and a broader range of
$B$
may reveal different behaviours (Zhao et al. Reference Zhao, Zhao, Zhou and Chong2024).
