Recent microfluidic experiments have explored the precise positioning of micron-sized particles in liquid environments via laser-induced thermoviscous flow. From micro-robotics to biology at the subcellular scale, this versatile technique has found a broad range of applications. Through the interplay between thermal expansion and thermal viscosity changes, the repeated scanning of the laser along a scan path results in fluid flow and hence net transport of particles, without physical channels. Building on previous work focusing on two-dimensional microfluidic settings, we present an analytical, theoretical model for the thermoviscous and thermal expansion flows and net transport induced by a translating heat spot in three-dimensional, unconfined fluid. We first numerically solve for the temperature field due to a translating heat source in the experimentally relevant limit. Then, in our flow model, the small, localised temperature increase causes local changes in the mass density, shear viscosity and bulk viscosity of the fluid. We derive analytically the instantaneous flow generated during one scan and compute the net transport of passive tracers due to a full scan, up to quadratic order in the thermal expansion and thermal shear viscosity coefficients. We further show that the flow and transport are independent of bulk viscosity. In the far field, while the leading-order instantaneous flow is typically a three-dimensional source or sink, the leading-order average velocity of tracers is instead a source dipole, whose strength depends on the relative magnitudes of the thermal expansion and thermal shear viscosity coefficients. Our quantitative results reveal the potential for future three-dimensional net transport and manipulation of particles at the microscale.
KeywordsCytoplasmic streamingMicrofluidicsStokes flowThermal expansionThermoviscous flowMathematics Subject Classification76-1076D0776N06EngineeringInterdisciplinary EngineeringEngineering and Physical Sciences Research Councilhttp://dx.doi.org/10.13039/501100000266Trinity College, University of Cambridgehttp://dx.doi.org/10.13039/501100000727publisher-imprint-nameSpringervolume-issue-count1issue-article-count6issue-toc-levels0issue-pricelist-year2025issue-copyright-holderSpringer Nature B.V.issue-copyright-year2025article-contains-esmNoarticle-numbering-styleContentOnlyarticle-registration-date-year2025article-registration-date-month3article-registration-date-day22article-toc-levels0toc-levels0volume-typeRegularjournal-productArchiveJournalnumbering-styleContentOnlyarticle-grants-typeOpenChoicemetadata-grantOpenAccessabstract-grantOpenAccessbodypdf-grantOpenAccessbodyhtml-grantOpenAccessbibliography-grantOpenAccessesm-grantOpenAccessonline-firstfalsepdf-file-referenceBodyRef/PDF/10665_2025_Article_10445.pdfpdf-typeTypesettarget-typeOnlinePDFissue-typeRegulararticle-typeOriginalPaperjournal-subject-primaryMathematicsjournal-subject-secondaryApplications of Mathematicsjournal-subject-secondaryMathematical Modeling and Industrial Mathematicsjournal-subject-secondaryComputational Mathematics and Numerical Analysisjournal-subject-secondaryMathematical and Computational Engineeringjournal-subject-secondaryTheoretical and Applied Mechanicsjournal-subject-collectionMathematics and Statisticsopen-accesstrueIntroduction
The controlled, precise manipulation of microparticles and biological cells in a liquid environment underpins applications across a wide range of disciplines and industries, from drug delivery and disease diagnostics to the construction of miniaturised devices and the study of single cells and molecules [1, 2–3]. Microparticle manipulation has been accomplished using a variety of physical mechanisms to achieve trapping, particle movement or assembly of structures [3]. Optical tweezers trap particles with optical forces associated with spatial gradients in light intensity, created by focused light [4, 5–6], while magnetic tweezers apply force to magnetic beads via a magnetic field [7], and acoustic tweezers use the acoustic radiation force induced by standing waves [8]. For optothermal micromanipulation, thermal gradients caused by laser heating result in thermophoresis, the migration of objects along a temperature gradient, providing control over microparticles [9]; active optothermophoretic manipulation involves particles that can establish their own temperature gradients [9].
With varied possibilities for driving mechanisms, fluid flow is instrumental in many techniques for micromanipulation. One of its advantages is that it need not rely on the optical or magnetic properties of the particles being manipulated; instead, recent studies have employed a vast range of flow field topologies in order to trap or manoeuvre particles hydrodynamically. For instance, time-averaged flow created by fluid oscillation, known as steady streaming, has enabled contact-free trapping of single cells in microeddies [8, 10]. Both driven by laser-induced temperature gradients, thermal convection and thermophoresis have been combined to manipulate microparticles [11], while Marangoni flows can be controlled optothermally [12] or via photoresponsive surfactants [13], to trap and guide the movement of particles at a gas-liquid interface, or to transport living cells [14]. Electroosmotic flow, which is actuated by electrodes that generate electric fields and hence flow of electrolyte solutions, has been used together with feedback control to independently steer and trap multiple particles at once [15]. In a microfluidic device with multiple channels, pressure-driven flow can be controlled via the channel flow rates; using this flow control, recent experiments have demonstrated precise manipulation of particles along user-defined paths, exploiting different flow modes [16, 17, 18, 19–20]. However, increasing the number of particles increases the number of degrees of freedom, necessitating a larger number of channels of the microfluidic device to control the particles independently [18].
Recent microfluidic experiments have investigated the use of thermoviscous flows for high-precision positioning of particles [21], trapping [22] and assembly [23]. In these studies, thermoviscous flow is driven by spatio-temporally varying heating of fluid, induced by a laser scanning along specified paths at kilohertz frequency. Although this technique, too, involves laser-induced heating, the physical mechanism for flow distinguishes it from the methods previously described. Here, the repeated scanning of the laser along a scan path results in fluid flow through the interplay between thermal expansion and thermal viscosity changes. This enables localised net transport of particles in the bulk fluid, without requiring direct laser exposure or physical channels [23, 24].
Highly versatile, these thermoviscous flows have also been applied inside living cells, and are known as focused-light-induced cytoplasmic streaming (FLUCS) in this context [25, 26–27]. This name contrasts FLUCS with naturally occurring cytoplasmic streaming: actively driven flows of the water-based, complex fluid (cytoplasm) inside cells, found in a wide variety of living organisms [28, 29]. For example, with thermoviscous flow perturbations, experiments have investigated how intracellular flows drive cell polarisation [25], and measured the rheology of cytoplasm [25] and the nucleus [26].
Previous theoretical work [24, 30, 31] has focused on the setup relevant to many of these microfluidic experiments [21, 22, 23, 24–25, 27, 30, 31], namely, with viscous fluid confined between parallel rigid plates. In Ref. [30], we considered a temperature profile with circular symmetry induced by the laser, which translates along a finite scan path, always in the same direction. We modelled the density and shear viscosity of the fluid as depending linearly on temperature change, a valid approximation for experimentally relevant temperature increases. The confinement provided by the parallel plates reduced the problem to two spatial dimensions, with a parabolic flow profile along the third dimension in the lubrication limit. We solved for the instantaneous flow field, and hence the average velocity of tracer particles throughout space due to repeated scanning. This revealed strong net transport near the scan path, together with inverse-square-law spatial decay, in agreement with experimental data [21]. Perhaps counter-intuitively, this net transport of tracers depends not on the magnitude, but instead on the rate of change of shear viscosity with temperature; hence, net thermoviscous flows are effective in fluids with vastly different viscosities [25, 30, 31]. Extending the method from one scan path to many, feedback algorithms have enabled the simultaneous manipulation of multiple particles, via net flow fields with complex topology [21, 23, 27], as shown in Fig. 1A and B.
Two-dimensional (A, B) and three-dimensional (C) thermoviscous and thermal expansion-driven net flows. A Trajectories of tracers in viscous fluid confined between parallel plates, induced by repeated scanning of six scan paths, in an experiment [23] (top) and according to analytical modelling of net thermoviscous flows [23, 30] (bottom). Scale bar: 15μm. B Control of 15 microparticles (white) with 8 scan paths (magenta) to form a humanoid figure in an experiment [23]. The net thermoviscous flow predicted by theory [23, 30] is shown in blue, with target positions in green. Scale bar: 15μm. A and B adapted from Ref. [23] and licensed under CC BY 4.0. C Theoretical trajectories of tracers in three-dimensional, unconfined fluid induced by scanning of a spherical heat spot (scan direction indicated by arrow), as derived in this work, due to thermoviscous effects (top) and due to purely thermal expansion-driven flows (bottom)
In general, micromanipulation can be conducted in two [7, 9, 11, 16, 17–18] or three [4, 7, 8–9, 19, 20, 32] dimensions. So far, thermoviscous flows have been used for two-dimensional micromanipulation experimentally [21, 22, 23–24], while existing theoretical studies consider essentially two-dimensional flow in the parallel-plate geometry [30, 31]. However, the possibility of three-dimensional thermoviscous and thermal expansion-driven net flows and micromanipulation, for example, with highly focused heating [24, 30], has yet to be explored. Unlocking more degrees of freedom could result in increased versatility of the technique for future experiments.
Geometry has a significant influence on the flows and transport induced by the focused light. For example, it can strongly modify the far-field behaviour, which, importantly, determines the impact of directing a particle towards its target on other particles, in microfluidic experiments [23]. Thus, to make progress towards three-dimensional micromanipulation, fundamental, quantitative understanding of three-dimensional flow is required. In this article, we study thermoviscous and thermal expansion-driven flows of three-dimensional, unconfined fluid.
As a consequence of the unconfined geometry, we must also include a qualitatively new physical ingredient, known as the bulk viscosity (or volume viscosity [33]) of the fluid, in our theory. This material property is relevant only for compressible flows. For the parallel-plate geometry, it could be safely neglected [30, 31]; however, without the parallel plates, bulk viscosity could play a role. From a continuum perspective, just like shear viscosity, bulk viscosity is a phenomenological coefficient that characterises a Newtonian fluid, and is independent of its state of motion [33]. While shear viscosity relates stress to linear deformation rate, bulk viscosity relates stress to volumetric deformation rate [33], i.e. dilatational-compressional motion of the fluid described by divergence of flow velocity [34]. At the molecular level, bulk viscosity reflects rotational and vibrational degrees of freedom in molecular motion, while shear viscosity is associated with translational motion [35, 36].
The bulk viscosity of a liquid may be determined experimentally via various techniques [34, 36, 37–38]. Motivated by recent microfluidic thermoviscous flow experiments [24, 25, 31], we may focus here on both glycerol-water mixtures and pure water. The authors of Ref. [37] demonstrated that the bulk viscosity of aqueous glycerol is comparable in magnitude with the shear viscosity, and that it depends on temperature. Furthermore, different groups have shown that the bulk viscosity of water decreases with temperature, at a rate comparable to that of shear viscosity, over the temperature range relevant to thermoviscous flow experiments [34, 37, 38]. Therefore, we cannot neglect the magnitude and rate of change with temperature of bulk viscosity in favour of those of shear viscosity, and here we retain both bulk viscosity and its temperature dependence in our theory.
In this paper, we solve systematically for the flow and net transport of particles that result from the scanning of a model spherical heat spot in unconfined viscous fluid. This article is organised as follows. We begin in Sect. 2 with the heat transport problem, solving numerically for the temperature field induced by a prescribed translating heat source in the experimentally relevant limit. Based on our results, we then choose a temperature profile that will serve as an input to our fluid flow model, representing the effect of the spatio-temporally varying heating. Next, in Sect. 3, we introduce our model for the compressible fluid flow, including the effects of thermal expansion, thermal shear viscosity changes and bulk viscosity with arbitrary temperature dependence. We derive an analytical expression for the flow induced during the translation of the heat spot up to quadratic order. Although our model includes the effect of bulk viscosity, we demonstrate that the flow is independent of it, instead depending on thermal expansion and thermal shear viscosity changes, just as for the parallel-plate geometry. We then compute in Sect. 4 the net transport of passive tracers induced by the flow, due to a full scan of the heat spot. In contrast with the instantaneous fluid flow, the leading-order net transport occurs at quadratic order in the thermal expansion coefficient and thermal shear viscosity coefficient. The direction of the net transport generated depends on the relative importance of these two effects, which can provide competing contributions: one from the interplay of thermal expansion and thermal shear viscosity changes (thermoviscous net transport, shown in the top panel of Fig. 1C), and the other purely driven by thermal expansion (Fig. 1C, bottom). Furthermore, the average velocity of tracers decays as a hydrodynamic source dipole in the far field, more strongly than the instantaneous flow during one scan of the heat spot. Finally, in Sect. 5, we compare the results of our axisymmetric model with our previous work for the parallel-plate setup, and propose future experiments to validate our theory, along with potential applications.
Temperature field induced by scanning heat source
In this section, we solve numerically for the temperature field induced by a prescribed heat source that moves along a scan path, in three-dimensional, unbounded space, in the limit relevant experimentally to FLUCS [25, 27]. Motivated by our results, we then introduce the temperature profile that we will prescribe in our model of the flow induced by the scanning heat source (Sect. 3).
Setup
Setup for our model of heat transport induced by a scanning laser. A spherical heat source of characteristic radius b translates at speed U in the z direction, along a scan path from z=-ℓ to z=ℓ along the z axis (cylindrical radial coordinate r), in unbounded, viscous fluid, causing a localised temperature perturbation
We illustrate the setup in Fig. 2, as a model for the standard FLUCS setup [25]. A prescribed heat source of characteristic radius b translates at constant speed U, along a scan path on the z axis from z=-ℓ to z=ℓ, in three-dimensional, unconfined fluid. The setup is axisymmetric about the z axis; we denote the cylindrical radial coordinate by r. We consider one scan of the heat source: the scan starts at time t=-t0 and finishes at time t=t0, so that the centre of the heat source has position (r=0,z=Ut) at time t, and the scan period is given by 2t0≡2ℓ/U. This heat source induces a localised temperature increase in the fluid, which in turn drives fluid flow (modelled in Sect. 3).
Governing equation for heat transport
We consider the heat transport problem that determines the temperature field T due to the prescribed heat source. We first make the simplifying assumption to neglect advection of heat by fluid flow: it may be verified a posteriori (Sect. 3.5) that the characteristic scale for the term ∂T/∂t is much larger than for the advective term u·∇T [30], where u is the fluid velocity field, under the same assumptions we employ to solve for flow; the key is that asymptotically, the flow velocity we will obtain in Sect. 3 is much smaller than the translation speed of the heat source. The temperature field therefore evolves primarily because of thermal diffusion and the prescribed heat source, independently of the flow driven by the temperature perturbation. Further, as a first approximation in this section on heat transport, we treat the material properties of the fluid as constant; however, when we consider the fluid flow driven by the temperature field in Sect. 3, we will allow the material properties to vary with temperature. Thus, the temperature field is governed by the forced heat equation asρ0cp∂T∂t=k∇2T+Φ,where ρ0, cp and k are the (constant) density, specific heat capacity and thermal conductivity of the fluid, respectively, and Φ is the source term, which provides the forcing. For the source term, we prescribe a translating, spherically symmetrical, Gaussian profile with time-dependent (non-negative) amplitude asΦ=Φ02πB(t)exp{-[r2+(z-Ut)2]/2b2},where Φ0 is the characteristic scale for the heat source term and the dimensionless amplitude function B(t) is given byB(t)=cosπt2t02for-t0≤t≤t0,0otherwise.Highly focused heating of the fluid has been suggested as a potential method of achieving this experimentally [24, 25]. Here, we consider only one scan, after which the heat source switches off completely, to investigate the temperature field (we will consider the effect of repeated scanning on net fluid flow in later sections). For convenience, we may write the temperature field T(r, z, t) of the fluid asT(r,z,t)=T0+ΔT(r,z,t),where T0 is a constant reference temperature and ΔT(r,z,t) is the temperature change of the fluid due to the localised heat source, which we assume to decay to zero at infinity.
Nondimensionalisation
For the temperature problem, we nondimensionalise length with b, time with b/U and temperature change with Φ0b2/k; we note that for the later sections on fluid flow (Sects. 3 and 4), we will use a different nondimensionalisation. In what follows, we use variable names to mean their dimensionless equivalents, to simplify notation. The dimensionless forced heat equation then becomesPescan∂ΔT∂t=∇2ΔT+Φ,where the dimensionless heat source term Φ is given byΦ=12πB(t)exp{-[r2+(z-t)2]/2},and we define the parameter Pescan to be the scanning Péclet number, given byPescan≡Ub(k/ρ0cp),where the denominator is the thermal diffusivity. Pescan is a dimensionless ratio that quantifies the relative importance of the scanning speed and thermal diffusion. It may be thought of as an unsteady Péclet number; importantly, the characteristic speed involved is the scanning speed, and not the fluid flow speed that would feature in the classical Péclet number.
Numerical simulation details
To obtain an estimate for Pescan, we substitute parameter values for water at 20∘C at atmospheric pressure [39], given by k=0.5861Wm-1K-1, ρ0=998.2kgm-3 and cp=4183Jkg-1K-1. Although the geometry is new in the present work, we use sample parameter values (dimensional) for the scan properties based on previous FLUCS and microfluidic experiments in confined geometries [21, 25], setting the heat source radius as b=4μm and the scan-path length as 2ℓ=11μm, so that the dimensionless half-scan period is given by t0=1.375. With these parameter values and selected dimensional scan periods, we solve the dimensionless heat equation (Eq. (5)) numerically, using finite element analysis in MATLAB R2024a. We present results for two different scanning Péclet numbers: first, Pescan=0.63, which corresponds to scan frequency 2k (i.e. dimensional scan period 2t0=0.5ms and speed U=0.022ms-1), typical for FLUCS experiments [21, 25, 27]; then Pescan=2, representing a faster scan within the experimental range.
Temperature profile for numerical simulations of forced heat equation for scanning Péclet number Pescan=0.63, during one scan of the heat source (-1.375≤t≤1.375). A Heat map showing spatial variation of temperature perturbation ΔT with cylindrical coordinates r and z, at selected times. The centre of the heat source is indicated in white, with the scan path in black. B Amplitude of temperature perturbation A(t) as a function of time. C Shape of temperature perturbation ΔT/A(t), i.e. temperature perturbation normalised by its peak, along the z axis, at selected times. The location of the peak temperature on the z axis is given by z=zpeak(t). Colours change from pink to yellow as time t increases. The dashed curve indicates the shape of the prescribed heat source (Eq. (8)), for comparison
First, for the slower scan of the heat source with scanning Péclet number Pescan=0.63 (parameter values in Sect. 2.4), we illustrate in Fig. 3 the temperature field induced during one scan period (-1.375≤t≤1.375), as obtained by numerical simulation of the forced heat equation. In Fig. 3A, we show snapshots of the temperature field as the heat source translates along the scan path. We observe that the region of higher temperature is highly localised and approximately spherically symmetrical. Like the prescribed heat source, this temperature perturbation translates in the positive z direction, first increasing in amplitude, peaking during the scan period and then decreasing. Any heat that remains at the end of the scan of the heat source diffuses away afterwards, so the temperature decays to ambient.
To further characterise the temperature field, we plot in Fig. 3B and Fig. 3C the amplitude and the shape of the temperature perturbation, respectively. Here, we define the amplitude A(t) (dimensionless) as the instantaneous maximum temperature perturbation at time t, which occurs at position z=zpeak(t) on the z axis, while the shape is the temperature perturbation field ΔT divided by this amplitude, given by ΔT/A(t). Over the course of a scan period, the amplitude (Fig. 3B) increases, peaks just after halfway through the scan period, and then decreases. Note also that the amplitude is not precisely zero at the end of the scan period, as it takes time for the heat remaining after the heat source has been switched off to diffuse away. The amplitude of the temperature perturbation therefore inherits broad qualitative features from that of the heat source, with differences in finer details.
In Fig. 3C, we plot the shape of the temperature perturbation along the z axis, as a function of displacement z-zpeak(t) from the location of the instantaneous maximum temperature, at selected times. For comparison, the black dashed curve is given by the formulaΔTA(t)=exp{-[z-zpeak(t)]2/2}.This is a Gaussian of the same radius as the prescribed heat source (centred at zpeak(t)). The key observation here is that strikingly, the shape is well approximated by a Gaussian throughout the scan. As time progresses, the shape exhibits a small increase in radius due to diffusive spreading. The characteristic distance by which the temperature perturbation spreads out can be determined via an intuitive scaling argument. In the (dimensional) scan period 2t0, heat diffuses by a characteristic distance of (k/ρ0cp)(2t0)≈8μm for our parameter values, which corresponds to a dimensionless distance of approximately 2, in good agreement with the degree of spreading seen in Fig. 3C.
Shape of temperature perturbation ΔT/A(t) along the z axis, at selected times, as obtained by numerical simulation of the forced heat equation for scanning Péclet number Pescan=2, during one scan of the heat source (-1.375≤t≤1.375), for comparison with result for Pescan=0.63 in Fig. 3C. The dashed curve shows the shape of the prescribed heat source (Eq. (8))
To explore the diffusive spreading of the Gaussian shape further, we now consider a higher scanning Péclet number (faster scanning). For Pescan=2, we illustrate the shape of the temperature perturbation ΔT/A(t) in Fig. 4. We see reduced spreading in this case. Indeed, in the limit of high Pescan, completely neglecting the Laplacian (diffusion) term on the right-hand side of the heat equation (Eq. (5)) allows us to integrate once with respect to time and hence write down the temperature field in this regime asΔT=12πPescanexp(-r2/2)∫-t0tB(t~)exp[-(z-t~)2/2]dt~.It may be verified that this essentially inherits the Gaussian shape and radius of the prescribed heat source (by inspection for the r dependence and by considering the far field for the z dependence).
Therefore, the initial Gaussian shape and radius originate from the heat source directly setting the evolution of the temperature field; during the scan, the shape remains Gaussian and the radius increases a little due to diffusion. Crucially, from a mathematical perspective, the unsteady term Pescan∂ΔT∂t in Eq. (5) must be retained. This is confirmed in the case of the two-dimensional (parallel-plate) geometry by experimental measurements showing that the temperature perturbation is highly localised [24, 25], an observation only consistent with the forced diffusion equation when the unsteady term is not discarded.
Model for temperature profile
We are now ready to propose a simplified temperature profile, as an input for our fluid flow model in the next section. To model the effect of the heating on the flow during one scan, we will prescribe a Gaussian temperature perturbation with time-varying amplitude, given dimensionally byΔT(r,z,t)=ΔT0A(t)exp{-[r2+(z-Ut)2]/2a2},where ΔT0 is the characteristic temperature change (a positive constant), A(t) is a dimensionless amplitude function (taking non-negative values), U is the speed of translation and a is the characteristic radius. Here, the amplitude function can be thought of as arbitrary or prescribed. This temperature profile captures the essential features of the solution to the heat transport problem we solved numerically, needed to explain thermal expansion-driven and thermoviscous fluid flows: the translation, shape and time-dependent amplitude. In particular, we have approximated the radius of the Gaussian as constant, neglecting the limited amount of diffusive spreading observed in the numerical simulations; in two dimensions, the thermoviscous net flows calculated theoretically under the same assumption have been shown to agree well with experimental results [23, 30].
In our modelling in Eq. (10), we need not select the same value for the characteristic radius of the temperature perturbation a as for that of the heat source b, although we may choose to as a modelling assumption. Similar applies for the amplitude function. We note also that we will later assume that the amplitude of the temperature perturbation is zero at the start and end of a scan, following Ref. [30]. This is a simplifying approximation, as we saw that the temperature decays to ambient over time at the end of a scan by diffusion, according to the heat equation. This approximation will allow us to treat every scan in repeated scanning as the same, as we assume that the temperature field resets to ambient at the end of every scan.
Finally, we recall that in the parallel-plate geometry, the temperature profile was also modelled as Gaussian [24, 30] based on experiments [24, 25], with cylindrical instead of spherical symmetry. Our analogous temperature profile for three-dimensional geometry will enable direct comparison of our results with the two-dimensional case.
Instantaneous flow during one scan
In this section, we introduce a three-dimensional model for the thermoviscous and thermal expansion-driven flows induced during one scan of a translating heat spot, in unbounded fluid. Here and in what follows, the term “heat spot” refers to the model temperature perturbation detailed in Sect. 2.6. With the instantaneous flow during one scan, we will compute in Sect. 4 the net transport of particles due to a full scan, which, by analogy with the two-dimensional case [30], we anticipate will allow trapping and manipulation of micron-sized particles in future experiments.
Setup
The setup is essentially as illustrated in Fig. 2. We now consider a localised, spherically symmetric temperature perturbation of characteristic radius a that translates along a scan path, given by the line segment from z=-ℓ to z=ℓ along the z axis, at speed U in unbounded, viscous fluid. During the scan, the centre of the heat spot is therefore at (r=0,z=Ut) at time t, for -t0≤t≤t0, where the scan period is 2t0≡2ℓ/U; the heat spot travels upwards. Again, the setup is axisymmetric about the z axis, with cylindrical radial coordinate denoted by r. We aim in this section to solve for the instantaneous fluid flow during one scan, driven by thermal expansion and thermal viscosity changes; we will build on this in Sect. 4 to understand the net transport of tracers due to repeated scanning of the heat spot, always upwards.
Governing equations
The translating localised temperature increase in the fluid modifies its material properties locally, driving fluid flow governed by the mass conservation and the momentum equations. Most of the explicit formulae for flow derived in this article correspond to a prescribed Gaussian temperature perturbation ΔT(r,z,t) with time-varying dimensionless amplitude A(t), introduced in Sect. 2 and given by Eq. (10). This model will allow us to understand the physics of thermoviscous and thermal expansion-driven flows. For small temperature changes, we use a standard linear relationship between the density of the fluid and the temperature, given byρ=ρ0(1-αΔT),where ρ0 is the density of the fluid at the reference temperature T0, and α is the thermal expansion coefficient (typically positive).
The mass conservation equation is given by∂ρ∂t+∇·(ρu)=0,where u is the velocity field; thus, the spatio-temporally varying density field gives rise to compressible fluid flow. The Cauchy momentum equation is given byρDuDt=∇·Π+ρg,where g is the gravitational acceleration and the stress tensor Π is given byΠ=-p1+κ(∇·u)1+2η12[∇u+(∇u)T]-131(∇·u),where p is the pressure field, 1 is the identity tensor, and η and κ are, respectively, the shear and bulk viscosities of the fluid. In this work, we consider a regime in which inertia and gravity are both negligible, relevant at small length scales; the dimensionless numbers associated with these conditions may be found by scaling arguments similar to those presented in Ref. [30]. For example, adapting to our geometry, the gravity-driven flow (natural convection) scales as ρ0ga2αΔT0/η0, whereas the thermal expansion-driven flow scales as αΔT0U. Thus, natural convection may be neglected if the dimensionless ratio ρ0ga2/η0U is small; we calculate its value as less than 0.01, for water at 20∘C at atmospheric pressure [39] and estimated parameter values of heat-spot radius a=4μm and speed U=0.022ms-1 [21, 25].
The shear viscosity η varies with temperature of the fluid, which we model (as in previous studies [24, 30, 31]) asη=η0(1-βΔT),where η0 is the shear viscosity of the fluid at the reference temperature T0, and β is the thermal shear viscosity coefficient (typically positive).
We recall that in the lubrication limit for the parallel-plate setup, the effect of bulk viscosity was shown to be unimportant [30, 31]. However, for aqueous glycerol, shear viscosity and bulk viscosity are similar in magnitude [34, 36, 37]; furthermore, bulk viscosity varies with temperature [34, 37, 38]. Thus, in our theory, we must retain the bulk viscosity. It would be possible to treat the temperature dependence of bulk viscosity in the same way as for density or shear viscosity for small temperature changes. That is, we could model the bulk viscosity as depending linearly on temperature, introducing a parameter to characterise the rate of variation. However, here, this turns out to be unnecessary and, as we show below, we can instead treat the bulk viscosity in full generality, without assuming a specific functional form for its dependence on temperature. We finally note that for a temperature perturbation that decays at infinity, the bulk viscosity, correspondingly, tends to a constant at infinity, just as the fluid density and shear viscosity do.
Nondimensionalisation
We nondimensionalise length with a, velocity with U, time with a/U, pressure with the viscous scale η0U/a, both shear viscosity and bulk viscosity with η0, density with ρ0 and temperature with ΔT0; note that this nondimensionalisation differs to that for the thermal problem in Sect. 2. In what follows, we use variable names to mean their dimensionless versions for notational convenience. We summarise the dimensionless governing equations. The momentum equation (Eqs. (13)–(14)) becomes-∇p+∇(κ∇·u)+η∇2u+13η∇(∇·u)+(∇η)·[∇u+(∇u)T]-23(∇η)(∇·u)=0.Here, we draw attention to the bulk viscosity term ∇(κ∇·u), which appears in the momentum equation as an exact gradient. Thus, a solution (u′,p′) to the governing equations with bulk viscosity set to zero induces a solution (u,p) to the full equations with arbitrary bulk viscosity, with the two pressure fields related via p′=p-κ∇·u.
After nondimensionalisation, the mass conservation equation (Eq. (12)) is still given by∂ρ∂t+∇·(ρu)=0.The density (Eq. (11)) and shear viscosity (Eq. (15)) of the fluid are now related to the temperature change byρ=1-αΔT,andη=1-βΔT,respectively, while the Gaussian model temperature perturbation (from Eq. (10)) becomesΔT(r,z,t)=A(t)exp{-[r2+(z-t)2]/2}.
Perturbation expansion
Due to our assumption of small temperature changes, the dimensionless coefficients α and β are now small parameters, representing the proportion by which the density and shear viscosity decrease in response to a temperature increase of ΔT0, respectively. We therefore solve this problem perturbatively, with the boundary conditions that the velocity and pressure are non-singular at the origin and decay at infinity (for temperature profiles that decay at infinity). We pose perturbation expansions for the velocity field and pressure field in the two small parameters, given byu=u0,0+αu1,0+βu0,1+α2u2,0+αβu1,1+β2u0,2+⋯,p=p0,0+αp1,0+βp0,1+α2p2,0+αβp1,1+β2p0,2+⋯,i.e. the velocity and pressure at order αmβn are given by um,n and pm,n, respectively. We note that had we chosen to model the bulk viscosity as depending linearly on temperature, we would have posed above a perturbation expansion in three small parameters instead of two, with the third parameter being a thermal bulk viscosity coefficient. In the following sections, we will solve for the flow at each order.
Solution at order <inline-formula id="IEq101"><alternatives><mml:math id="IEq101_Math"><mml:msup><mml:mi>β</mml:mi><mml:mi>n</mml:mi></mml:msup></mml:math><tex-math id="IEq101_TeX"><![CDATA[\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\beta ^{n}$$\end{document}]]></tex-math><inline-graphic specific-use="web" mime-subtype="GIF" xlink:href="10665_2025_10445_Article_IEq101.gif"/></alternatives></inline-formula>
First, as was shown for the parallel-plate geometry [30], we claim that in Eqs. (21) and (22), we may set both the velocity and pressure to be zero at orders βn for all n, with the physical interpretation that thermal expansion is essential for the fluid flow. To see this, we observe that if we set the thermal expansion coefficient α to be zero, then the terms at order αmβn with m≥1 in Eqs. (21) and (22) vanish, while the mass conservation equation (Eq. (17)) becomes∇·u=0.We observe that zero flow, u=0 and p=0, solves Eq. (23) together with Eq. (16). By expanding this solution in the parameter β, we then see that this corresponds to u0,n=0 and p0,n=0, as claimed. The perturbation expansions then simplify tou=αu1,0+α2u2,0+αβu1,1+⋯,p=αp1,0+α2p2,0+αβp1,1+⋯.The structure of these expansions reflects the physics of the flows. Every term includes the thermal expansion coefficient α, which mathematically captures the fact that thermal expansion drives the flow, via forcing in the mass conservation equation. This flow is an automatic consequence of the spatio-temporally varying temperature field [40]; we emphasise that it is independent of gravity, thus distinguishing it from buoyancy-driven flows.
We are now in a position to revisit the assumption to neglect advection of heat made in Sect. 2 for the thermal problem. Since the instantaneous flow occurs at order α, the ratio of the advection term u·∇T to the rate of change of temperature at a fixed position due to scanning ∂T/∂t is an order-α quantity; thus, we confirm that we may neglect advection of heat by fluid flow in favour of the unsteady term in the leading-order thermal problem.
Velocity field associated with the time variation of the heat-spot amplitude
We may decompose the fluid flow induced by the scanning heat spot into two contributions, in a precise manner that we explain mathematically below: one associated with the time variation of the heat-spot amplitude, and the other related to the translation of the heat spot. In this section, we consider the former and derive the contribution to the velocity field that captures the switching-on and switching-off of a spherically symmetric heat spot.
Decomposition of velocity field
First, we introduce spherical polar coordinates (R,θ,ϕ) with origin at the centre of the translating heat spot, i.e. at (r=0,z=t), so that the spherical radial coordinate is given byR=r2+(z-t)2.We assume that the temperature perturbation during one scan of the heat spot has the formΔT(R,t)=A(t)Θ(R),i.e. an arbitrary, time-dependent amplitude A(t), which is zero at the ends of the scan path, multiplied by a shape function Θ(R) that translates in the z direction and decays at infinity. The spherical symmetry of the (instantaneous) heat-spot shape is mathematically convenient for this first model in three dimensions. In this case, using Eqs. (18), (26) and (27), the mass conservation equation (Eq. (17)) becomes-αA′(t)Θ(R)+A(t)Θ′(R)∂R∂t+∇·(ρu)=0.The two forcing terms reflect the two reasons that the temperature field at a given position varies with time: because the amplitude is time-dependent (first term) and because the shape of the temperature perturbation translates in space (second term). Since Eq. (28) is linear in velocity, we may decompose the solution into two contributions (as in Ref. [30]), to account for the two forcing terms. We thus write the velocity field asu=u(S)+u(T).Here, we introduce the switching-on velocity field u(S), associated with the time variation of the heat-spot amplitude, as a solution to the equation given by-αA′(t)Θ(R)+∇·(ρu(S))=0,i.e. treating the forcing term containing A′(t), while the velocity field u(T) associated with translation of the heat spot satisfies-αA(t)Θ′(R)∂R∂t+∇·(ρu(T))=0,i.e. accounting for the forcing term containing Θ′(R)∂R∂t.
Mathematical derivation of switching-on flow
We treat in this section the switching-on velocity field, associated with the time variation of the heat-spot amplitude; we will return to the translational contribution in Sect. 3.7. The switching-on flow u(S) satisfies Eq. (30), a version of the mass conservation equation but forced by a regularised source, instantaneously centred at the location of the heat spot and with time-varying amplitude. We assume that this switching-on velocity field also satisfies the momentum equation with a corresponding switching-on pressure field, is non-singular at the origin and decays at infinity (provided the temperature perturbation decays sufficiently fast at infinity).
We now pose a spherically symmetric ansatz for the switching-on velocity field, given byu(S)=u(S)(R,t)eR,where eR is the radial unit vector from the centre of the heat spot. The mass conservation equation (Eq. (30)) in spherical polar coordinates then becomes-αA′(t)Θ(R)+1R2∂∂R(R2ρu(S))=0.
Flow result and physical interpretation
Integrating the defining mass conservation equation for u(S) (Eq. (33)) and applying the boundary conditions (Sect. 3.6.2), we obtain the switching-on flow (Eq. (32)), for an arbitrary, spherically symmetric temperature perturbation of the form in Eq. (27), asu(S)(R,t)=αA′(t)ρ(R,t)R2∫0RR~2Θ(R~)dR~.We observe that the perturbation expansion for this switching-on velocity u(S) features only powers of the thermal expansion coefficient α (from Taylor expanding the factor of ρ-1), but not the other dimensionless parameter β, which characterises thermal shear viscosity changes. It is also independent of the bulk viscosity κ. As in two dimensions [30], this switching-on contribution can be derived from only mass conservation and spherical symmetry of the prescribed temperature perturbation, independent of the momentum equation and hence independent of viscosity; it is a kinematic phenomenon that relies solely on thermal expansion. The flow u(S) also fully accounts for any appearances of the rate of change of heat-spot amplitude, A′(t), in the solution for the full flow u due to the spherically symmetric heat spot. Physically, when the heat-spot amplitude is increasing [A′(t)>00$$\end{document}]]>], the switching-on flow is an instantaneous regularised source, as the fluid expands radially outwards. Conversely, when the heat-spot amplitude is decreasing [A′(t)<0], the switching-on flow is instead an instantaneous sink, reflecting the contraction of the fluid as temperature decreases locally.
Solution at order <inline-formula id="IEq136"><alternatives><mml:math id="IEq136_Math"><mml:mi>α</mml:mi></mml:math><tex-math id="IEq136_TeX"><![CDATA[\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\alpha $$\end{document}]]></tex-math><inline-graphic specific-use="web" mime-subtype="GIF" xlink:href="10665_2025_10445_Article_IEq136.gif"/></alternatives></inline-formula>
We now proceed to solve in this section for the leading-order instantaneous flow during a scan, which occurs at order α and is purely driven by thermal expansion. We will show that this consists of a contribution due to the switching-on of the heat spot and another related to the translation of the heat spot.
General solution
With Eqs. (18) and (24) for the density and velocity, respectively, we expand Eq. (17) to obtain that at order α, the statement of mass conservation is given by-∂ΔT∂t+∇·u1,0=0,while the momentum equation (Eq. (16)) is given by-∇p1,0+∇κ+13∇·u1,0+∇2u1,0=0.Building on the decomposition of the flow (Eq. (29)) introduced in Sect. 3.6.1, we pose an ansatz for the velocity field at order α given byu1,0=∂∂t[A(t)u1,0(S)]≡A′(t)u1,0(S)+A(t)u1,0(T),for -t0≤t≤t0, where u1,0(S) is the velocity field associated with the time variation of the heat-spot amplitude and u1,0(T) is the velocity field associated with the translation of the heat spot. Its structure is inherited from the forcing term in the mass conservation equation (Eq. (35)), which takes the form of a time-derivative of a function proportional to the amplitude A(t). Here, we define the switching-on velocity field u1,0(S) at order α via a perturbation expansion of the switching-on velocity field u(S) (Eq. (34)) derived in Sect. 3.6, given byu(S)=αA′(t)R2∫0RR~2Θ(R~)dR~+O(α2)eR≡αA′(t)u1,0(S)+O(α2).It may be verified that the ansatz for the flow at order α satisfies the mass conservation and momentum equations, with the pressure field (found by taking the divergence of the momentum equation, i.e. Eq. (36), and combining with mass conservation) given byp1,0=κ+43∂ΔT∂t.We emphasise that we did not need to assume that the bulk viscosity κ is a constant; instead, it depends on space and time via the temperature field. While the pressure depends on bulk viscosity, we note that the velocity field at order α is independent of bulk viscosity. We also note that the velocity field at order α during a scan is the time-derivative of a function proportional to the heat-spot amplitude. Consequently, when integrated over a scan, this will give rise to zero net displacement of material points at order α, just as for the parallel-plate setup in Ref. [30].
Flow result for Gaussian temperature profile and physical interpretation
To illustrate our solution for the flow at order α during a scan given by Eqs. (37)–(38), we now prescribe a Gaussian temperature profile for the spherical heat spot, motivated in Sect. 2 and given byΔT(R,t)=A(t)exp(-R2/2).For this Gaussian temperature profile, the switching-on velocity field at order α, from Eq. (38), is given byu1,0(S)=1R2∫0RR~2exp(-R~2/2)dR~eR≡πerf(R/2)2R2-exp(-R2/2)ReR.This purely radial flow, shown in Fig. 5, is a hydrodynamic source in the far field, decaying as 1/R2 (since the fluid is three-dimensional and unbounded). As for the parallel-plate case in Ref. [30], this arises from local volume increase of the fluid due to heating.
Streamlines of the flow u1,0(S) at order α associated with the switching-on of the spherical heat spot, with centre at (r=0,z=t) in dimensionless coordinates (nondimensionalisation described in Sect. 3.3). Colour shows the magnitude |u1,0(S)|
The velocity field associated with the translation of the heat spot, from Eqs. (37) and (41), and shown in Fig. 6, is given byu1,0(T)≡∂∂tu1,0(S)=cosθ2πerf(R/2)R3-1+2R2exp(-R2/2)eR+sinθπerf(R/2)2R3-exp(-R2/2)R2eθ,where eθ is the spherical basis vector corresponding to the polar angle θ. The streamlines are qualitatively the same as those for the parallel-plate setup [30]. In the far field, the flow is a source dipole in the far field, decaying as 1/R3. The source at the front corresponds to the arrival of the heat spot, while the sink at the back results from cooling due to the departure of the heat spot, just as for the parallel-plate geometry [24, 30].
Streamlines of the instantaneous flow u1,0(T) at order α associated with the translation of the spherical heat spot, during one scan. The flow is axisymmetric about the z axis. Colour indicates the magnitude |u1,0(T)|
In terms of dimensional variables, the leading-order flow during one scan scales as αΔT0U; it is proportional to the thermal expansion coefficient, the peak temperature change and the velocity of translation of the heat spot. Hence, the leading-order instantaneous flow scales linearly with the scan frequency. Furthermore, changing the sign of the thermal expansion coefficient reverses this flow; for example, water has a positive value of α at room temperature, but negative when below 4∘C.
Note that for all plots of instantaneous fluid flow during one scan in Sect. 3, the vertical axis label is z-t, reflecting that the origin of our spherical polar coordinate system is at the centre of the heat spot (r=0,z=t). Thus, the plots show the flow as viewed when moving at the same velocity as the heat spot.
For water and glycerol (common liquids used in experiments), the thermal expansion coefficient α is much smaller than the thermal shear viscosity coefficient β, so we may expect effects at order αβ to be more visible in experiments using these fluids than those at order α2, as was the case for the parallel-plate setup [30]. We therefore now consider the solution for instantaneous thermoviscous flow during one scan at order αβ, and will next return to the purely thermal expansion-driven effect at order α2.
Mathematical derivation
At order αβ, mass conservation (Eq. (17)) is given by∇·u1,1=0,i.e. the flow at order αβ is incompressible. Here, we have substituted Eqs. (18) (density) and (24) (velocity perturbation expansion), and used the fact that there is no flow at order β (Sect. 3.5).
Expanding Eq. (16), the momentum equation at order αβ is given by-∇p1,1+∇(κ∇·u1,1)+∇2u1,1-ΔT∇2u1,0+13∇(∇·u1,1)-13ΔT∇(∇·u1,0)-(∇(ΔT))·[∇u1,0+(∇u1,0)T]+23(∇(ΔT))(∇·u1,0)=0.Using mass conservation at orders α (Eq. (35)) and αβ (Eq. (43)) to simplify this, we find-∇p1,1+∇2u1,1-ΔT∇2u1,0-13ΔT∇∂ΔT∂t-(∇(ΔT))·[∇u1,0+(∇u1,0)T]+23∇(ΔT)∂ΔT∂t=0.We observe that due to incompressibility of the flow at order αβ (Eq. (43)), the bulk viscosity κ does not feature in this equation, which takes the form of the forced Stokes equation; the solution for the velocity field u1,1 will correspondingly also be independent of bulk viscosity.
We now use a streamfunction approach to solve for the flow. The velocity field at order αβ is divergence-free, so we may write it in terms of a Stokes streamfunction Ψ asu1,1=∇×-Ψ(R,θ,t)Rsinθeϕ,where eϕ is the spherical basis vector corresponding to the azimuthal angle ϕ. It is a classical result [33] that∇2(∇×u1,1)=E2(E2Ψ)Rsinθeϕ,where the operator E2 is given byE2Ψ≡∂2Ψ∂R2+sinθR2∂∂θ1sinθ∂Ψ∂θ.Using this, the ϕ component of the curl of the momentum equation (Eq. (45)) in spherical polar coordinates becomesE2(E2Ψ)Rsinθ=A(t)2sinθ×4Rexp(-R2)+6Rexp(-R2)-32πR2exp(-R2/2)erf(R/2).Again, as with the parallel-plate setup [30], we observe that the forcing at order αβ does not contain A′(t); only the flow fields at orders αn (solely driven by thermal expansion) depend on the rate of change of the heat-spot amplitude, through the switching-on contribution.
To solve this, we choose the ansatzΨ=A(t)2f(R)sin2θ,which is similar to that for the incompressible Stokes flow past a rigid sphere [33]. By reducing the problem to two second-order ordinary differential equations in R (one for each application of the operator E2), this allows us to solve for the axisymmetric streamfunction Ψ asΨ=A(t)2sin2θ{π20Rerf(R)+π52R3+22πR-2π5Rexp(-R2/2)erf(R/2)+π20R4+π4R2[erf(R/2)2-1]+110R2+310exp(-R2)}.Applying Eq. (46) then gives the corresponding velocity field at order αβ, which satisfies the boundary conditions at the origin and at infinity.
Flow result and physical interpretation
The velocity field u1,1 at order αβ, derived above, is given in spherical polar coordinates byu1,1=A(t)2{eRcosθ{-π10R3erf(R)+-2π5R-42π5R+22π5R3exp(-R2/2)erf(R/2)+π10R2+π2[1-erf(R/2)2]-15+35R2exp(-R2)}+eθsinθ{-π20R3erf(R)+22π5R+32π5R+2π5R3exp(-R2/2)erf(R/2)+π5R2+π2[erf(R/2)2-1]+25-310R2exp(-R2)}},for -t0≤t≤t0. We plot the streamlines of this flow in Fig. 7 and note that it qualitatively matches that for the parallel-plate setup at the same order [30].
Streamlines of the flow u1,1 at order αβ due to a spherical heat spot, translating in the z direction. Colour shows the magnitude of the scaled flow |u1,1/A(t)2|
In the far field, the flow at order αβ is given byu1,1∼-A(t)2π20R3(2cosθeR+sinθeθ).This is a three-dimensional source dipole, decaying as 1/R3, with a sink at the front and a source at the back, the same physics as for the parallel-plate geometry [30]. Strikingly, the flow on the z axis is downwards, i.e. in the opposite direction to heat-spot translation. To explain the flow u1,1 at order αβ physically, we consider it as a modification, due to thermal shear viscosity changes, of the flow contribution A(t)u1,0(T) at order α associated with translation of the heat spot. Specifically, we may view the near field of u1,1 as localised amplification of the front-to-back flow near the heat spot of A(t)u1,0(T), due to reduced shear viscosity locally, while the far-field source dipole of u1,1 enforces incompressibility at this order.
Dimensionally, the thermoviscous flow at order αβ during one scan scales with αβΔT02U. Like the leading-order flow at order α, the order-αβ contribution is proportional to the velocity of the heat spot, and would be reversed by a change in sign of the thermal expansion coefficient. However, in contrast with the leading order, the flow at order αβ scales quadratically with peak temperature change.
We will see in Sect. 4 that this instantaneous thermoviscous flow at order αβ during a scan will contribute to the leading-order time-averaged velocity of tracers.
Solution at order <inline-formula id="IEq218"><alternatives><mml:math id="IEq218_Math"><mml:msup><mml:mi>α</mml:mi><mml:mn>2</mml:mn></mml:msup></mml:math><tex-math id="IEq218_TeX"><![CDATA[\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\alpha ^2$$\end{document}]]></tex-math><inline-graphic specific-use="web" mime-subtype="GIF" xlink:href="10665_2025_10445_Article_IEq218.gif"/></alternatives></inline-formula>
To complete the analysis of the instantaneous flow correct to quadratic order, we now consider the flow at order α2. This flow does not depend on thermal shear viscosity changes, instead relying only on thermal expansion; it is thus fundamentally different to the order-αβ thermoviscous flow contribution and would exist even if the shear viscosity were constant with respect to temperature.
Mathematical derivation
Expanding Eq. (17) (again using Eqs. (18) and (24)), we find that the mass conservation equation at order α2 is given by∇·u2,0-∇·(ΔTu1,0)=0.In contrast with order αβ, the flow at order α2 is not incompressible. However, to solve for the flow, it is useful to define an incompressible modified velocity field v asv≡u2,0-ΔTu1,0.From Eq. (16), at order α2, the momentum equation is given by-∇p2,0+∇(κ∇·u2,0)+∇2u2,0+13∇(∇·u2,0)=0.Using mass conservation (Eq. (54)), we simplify this to-∇p2,0+∇κ+13∇·(ΔTu1,0)+∇2u2,0=0.Taking the curl, we obtain∇2(∇×u2,0)=0,thus eliminating the gradient terms. In particular, the bulk viscosity term vanishes again, so that the solution for the flow u2,0 at order α2 will be independent of bulk viscosity. Rewriting this in terms of the modified velocity field (Eq. (55)) gives∇2[∇×(v+ΔTu1,0)]=0.The solution method now closely follows that at order αβ: since the modified velocity field is divergence-free, we may write it in terms of a Stokes streamfunction Φ asv=∇×-Φ(R,θ,t)Rsinθeϕ.Similarly to order αβ, we obtain the relation∇2(∇×v)=E2(E2Φ)Rsinθeϕ,where the operator E2 is given by Eq. (48). We then write Eq. (59) in terms of the streamfunction asE2(E2Φ)Rsinθeϕ=-∇2[∇×(ΔTu1,0)].The right-hand side of Eq. (62) is given explicitly by-∇2[∇×(ΔTu1,0)]=A(t)2sinθ[π21+1R2exp(-R2/2)erf(R/2)-4R+1Rexp(-R2/2)]eϕ.This does not contain the rate of change of the heat-spot amplitude A′(t), since the contribution associated with the switching-on of the heat spot involves the curl of a function of R parallel to the radial direction (ΔTu1,0(S)); this curl is zero by symmetry. As a result, we will see that the modified velocity field v does not contain A′(t). Any contributions to the velocity field u2,0 at order α2 involving the rate of change of heat-spot amplitude therefore originate from the term ΔTu1,0 in Eq. (55), specifically from ΔTu1,0(S).
Now, the streamfunction equation (Eq. (62)) becomesE2(E2Φ)=A(t)2sin2θπ2R+1Rexp(-R2/2)erf(R/2)-4R2+1exp(-R2/2).Following the same method as at order αβ, we choose the ansatzΦ=A(t)2g(R)sin2θ,and find that the function g(R) is given byg(R)=112(-πerf(R)R-22πRexp(-R2/2)erf(R/2)+22πexp(-R2/2)erf(R/2)R-2exp(-R2)+πR2-πR2erf(R/2)2).The Stokes streamfunction for the modified velocity field is therefore given by Eqs. (65)–(66). We can substitute this into Eq. (60) to find the modified velocity field, and hence obtain the flow u2,0 at order α2 from Eq. (55).
Flow result and physical interpretation
Streamlines of the radial flow u2,0(S) at order α2 associated with the switching-on of a spherical heat spot, with centre at (r=0,z=t). Colour shows the magnitude |u2,0(S)|
Streamlines of the axisymmetric flow u2,0(T) at order α2 due to a translating spherical heat spot during one scan. Colour indicates magnitude |u2,0(T)|
As detailed above, we obtain the solution for the velocity field at order α2 asu2,0≡A′(t)A(t)u2,0(S)+A(t)2u2,0(T),where the switching-on contribution is given byu2,0(S)=2πexp(-R2/2)erf(R/2)2R2-exp(-R2)ReR,and the translation contribution is given byu2,0(T)=cosθeR(πerf(R)6R3+2πexp(-R2/2)erf(R/2)3R+22πexp(-R2/2)erf(R/2)3R3-exp(-R2)-5exp(-R2)3R2+16π[erf(R/2)2-1])+sinθeθ(πerf(R)12R3-2πexp(-R2/2)erf(R/2)3R+2πexp(-R2/2)erf(R/2)3R3-5exp(-R2)6R2-16π[erf(R/2)2-1]),for -t0≤t≤t0 (i.e. during one scan). We illustrate the radial flow associated with the switching-on of the heat spot u2,0(S) in Fig. 8, and the dipolar flow associated with the translation of the heat spot u2,0(T) in Fig. 9. These two flows are reminiscent of and may be viewed as reinforcing their counterparts at order α; the flow at order α2 is the first correction to the flow at order α due to the fact that the density is lowered by heating below its reference value ρ0 for positive thermal expansion coefficient (α>00$$\end{document}]]>), in the flux term ∇·(ρu) in mass conservation (Eq. (12)).
The far-field behaviour of the flow u2,0 at order α2 is given byu2,0∼A(t)2π12R3(2cosθeR+sinθeθ).This is opposite in sign to the far field of order αβ, as for the parallel-plate setup [30]. However, deviating from the two-dimensional case [30], the numerical prefactor has a different magnitude compared with that at order αβ for the axisymmetric setup without confinement; although much of the physics is shared between the different geometries, our calculations demonstrate explicitly how the details of the flow are specific to each geometry.
In dimensional terms, the flow during a scan at order α2 scales as α2ΔT02U; in contrast with order αβ, changing the sign of the thermal expansion coefficient does not affect the contribution at order α2.
Mathematically, for a general fluid, we will show that this purely thermal expansion-driven flow at order α2 gives rise to a contribution to the leading-order average velocity of tracers over a scan, typically in the opposite direction to the thermoviscous contribution at order αβ.
Net transport of tracers
We introduced in Sect. 3 our model for the thermoviscous and thermal expansion-driven flows induced by a translating spherical heat spot in three-dimensional, unbounded fluid, solving for the instantaneous flow during one scan of the heat spot, up to quadratic order in the thermal expansion coefficient α and thermal shear viscosity coefficient β. In two-dimensional experiments, the relevant physical quantity observed is the net displacement of tracer beads due to repeated scanning of the laser, instead of the instantaneous fluid flow during one scan; furthermore, the theoretical average velocity of tracers from hydrodynamic modelling has been used quantitatively to design scan patterns to manipulate particles [23, 27, 30, 41]. We therefore anticipate that future three-dimensional experiments, like their two-dimensional counterparts, will measure and exploit the time-averaged trajectories of tracers; consequently, predicting the net displacement of tracers due to our three-dimensional, unconfined fluid flow over the course of one scan is a key step towards applications in micromanipulation. We thus now examine the kinematics of tracer particles in this flow and compute the net displacement of tracers due to a full scan of the heat spot from z=-ℓ to z=ℓ.
Trajectory of tracer
We solve for the net displacement of a material point due to one full scan of the heat spot along the scan path from z=-ℓ to z=ℓ. Consider a material point that has initial position X0 at time t=-t0. Following the method in Ref. [30], we write its position vector relative to the origin at time t as X(t). In the absence of noise, this obeys an ordinary differential equation given bydXdt=u(X(t),t),for -t0≤t≤t0 (i.e. during one scan). The equivalent integral equation is given byX(t)-X0=∫-t0tu(X(t~),t~)dt~.
Perturbation expansion
As in Ref. [30], we pose a perturbation expansion for the displacement vector ΔX(t)≡X(t)-X0 asΔX(t)=αΔX1,0(t)+α2ΔX2,0(t)+αβΔX1,1(t)+h.o.t.,where ΔXm,n(t) is the order-αmβn displacement of the material point at time t from the position X0 at t=-t0. Then, expanding Eq. (72) yieldsαΔX1,0(t)+α2ΔX2,0(t)+αβΔX1,1(t)+h.o.t.=α∫-t0tu1,0(X0,t~)dt~+α2∫-t0t[u2,0(X0,t~)+ΔX1,0(t~)·∇u1,0(X0,t~)]dt~+αβ∫-t0tu1,1(X0,t~)dt~.
Zero net displacement at order <inline-formula id="IEq285"><alternatives><mml:math id="IEq285_Math"><mml:mi>α</mml:mi></mml:math><tex-math id="IEq285_TeX"><![CDATA[\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\alpha $$\end{document}]]></tex-math><inline-graphic specific-use="web" mime-subtype="GIF" xlink:href="10665_2025_10445_Article_IEq285.gif"/></alternatives></inline-formula>
We show here that the net displacement of a tracer after one full scan of the heat spot varies not linearly with the temperature change, but instead (at least) quadratically. From the perturbation expansion in Eq. (74), the displacement ΔX1,0(t) at order α (i.e. leading order) of a material point at time t is given byΔX1,0(t)=∫-t0tu1,0(X0,t~)dt~.Using Eq. (37) for u1,0 and the Fundamental Theorem of Calculus, for a general heat spot, the expression for the order-α displacement in Eq. (75) becomesΔX1,0(t)=∫-t0t∂∂t~[A(t~)u1,0(S)(X0,t~)]dt~=A(t)u1,0(S)(X0,t)-A(-t0)u1,0(S)(X0,-t0).For a scan path of finite length, the heat-spot amplitude is zero at the ends of the scan path, so this simplifies toΔX1,0(t)=A(t)u1,0(S)(X0,t).Furthermore, the net displacement ΔX1,0(t0) (due to a full scan) at order α, of any material point, is given byΔX1,0(t0)=0.Thus, importantly, the leading-order net displacement of a tracer occurs not at linear order, but instead at quadratic order in the dimensionless parameters α and β; net transport is hence quadratic in the temperature perturbation.
Net displacement at order <inline-formula id="IEq294"><alternatives><mml:math id="IEq294_Math"><mml:mrow><mml:mi>α</mml:mi><mml:mi>β</mml:mi></mml:mrow></mml:math><tex-math id="IEq294_TeX"><![CDATA[\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\alpha \beta $$\end{document}]]></tex-math><inline-graphic specific-use="web" mime-subtype="GIF" xlink:href="10665_2025_10445_Article_IEq294.gif"/></alternatives></inline-formula> and order <inline-formula id="IEq295"><alternatives><mml:math id="IEq295_Math"><mml:msup><mml:mi>α</mml:mi><mml:mn>2</mml:mn></mml:msup></mml:math><tex-math id="IEq295_TeX"><![CDATA[\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\alpha ^2$$\end{document}]]></tex-math><inline-graphic specific-use="web" mime-subtype="GIF" xlink:href="10665_2025_10445_Article_IEq295.gif"/></alternatives></inline-formula>
Since we showed above that the net displacement at order α is precisely zero, the perturbation expansion for the net displacement of the material point (Eq. (73)) can be rewritten asΔX(t0)=α2ΔX2,0(t0)+αβΔX1,1(t0)+h.o.t.Here, by Eq. (74), the thermoviscous net displacement ΔX1,1(t0) at order αβ of a material point with initial position X0 is given byΔX1,1(t0)=∫-t0t0u1,1(X0,t)dt,and, similarly, the purely thermal expansion-driven net displacement ΔX2,0(t0) at order α2 is given byΔX2,0(t0)=∫-t0t0[u2,0(X0,t)+ΔX1,0(t)·∇u1,0(X0,t)]dt.
Net displacement ΔX1,1(t0) of tracers at order αβ due to one full scan of a spherical heat spot, with sinusoidal amplitude function (Eq. (82)), and scan path from z=-1.375 to z=1.375 along the z axis, with dimensionless scan-path length 2t0=2.75. Colour indicates the magnitude |ΔX1,1(t0)|
Net displacement ΔX2,0(t0) of tracers at order α2 due to one scan of a spherical heat spot from z=-1.375 to z=1.375. Colour shows the magnitude |ΔX2,0(t0)|
We plot the net displacement of tracers at order αβ (Eq. (80)) in Fig. 10 and at order α2 (Eq. (81)) in Fig. 11. Here, to facilitate direct comparison with the results of Ref. [30], we choose a sinusoidal heat-spot amplitude function, given byA(t)=cos2πt2t0,for -t0≤t≤t0, and we set the dimensionless scan-path length to be 2t0=2.75, so that the scan path is from z=-1.375 to z=1.375 along the z axis. We note that the vertical axis label for these plots showing net transport of tracers due to a full scan is z, instead of z-t used for the plots of instantaneous flow in Sect. 3; that is, the net transport is shown here in the laboratory frame, as opposed to the co-moving frame.
Net transport near the scan path
We now turn our attention to physical interpretation of our results on net transport (Eqs. (80) and (81)), illustrated in Figs. 10 and 11, beginning with the behaviour close to the scan path. Near the scan path, the upward translation of a heat spot thus produces two contributions to net transport of tracers, which typically have opposite directions, and scale dimensionally with αβΔT02 and α2ΔT02. First, assuming positive thermal expansion coefficient α and thermal shear viscosity coefficient β, the net transport at order αβ near the scan path is in the opposite direction to heat-spot translation (Fig. 10); this is an extension of the classic net thermoviscous flow of fluid confined between parallel plates and arises from the combination of thermal expansion and temperature-dependent shear viscosity [21, 24, 30, 31]. The physical mechanism for this thermoviscous net flow is inherited via the time-averaging in Eq. (80) from the instantaneous flow at order αβ, explained in Sect. 3.8.2. A negative thermal expansion coefficient α (e.g. for water below 4∘C) would reverse the direction of the contribution at order αβ to net displacement of tracers, by changing the sign of αβ in the perturbation expansion in Eq. (79).
Secondly, thermal expansion associated with the scanning heat spot can, by itself, give rise to nonzero net transport at order α2, via a fundamentally different physical mechanism. This is independent of thermal shear viscosity changes and is inherited from the flow at order α2 during one scan (Sect. 3.9). In principle, achieving net transport therefore does not require temperature-dependent viscosity, only thermal expansion, though in practice, its significance would vary based on the material parameters of the liquid used. This contribution is characterised by net transport always in the same direction as the heat-spot translation (upwards), whether the thermal expansion coefficient α is positive or negative. For aqueous glycerol, we expect the thermoviscous net transport at order αβ to dominate over the purely thermal expansion-driven contribution at order α2. However, if a liquid has sufficiently large thermal expansion coefficient relative to the thermal shear viscosity coefficient, then the net transport at order α2 may dominate instead, with opposite direction. In all cases, there is no contribution from the bulk viscosity.
Far-field behaviour
We finally examine the far-field behaviour of the thermoviscous and thermal expansion-driven net transport. In the far field, the net displacement at order αβ is given byΔX1,1(t0)∼-π20∫-t0t0A(t)2dt1R3(2cosθeR+sinθeθ)t=-t0,which is inherited from the source dipole in the far field of the flow u1,1 via Eqs. (80), (52) and (53), with physical interpretation given in Sect. 3.8.2. We note that evaluating at t=-t0 corresponds to spherical polar coordinates with origin at (r=0,z=-t0) (in terms of cylindrical coordinates), with the tracer at its initial position. Similarly, the far-field net displacement at order α2 is given byΔX2,0(t0)∼π12∫-t0t0A(t)2dt1R3(2cosθeR+sinθeθ)t=-t0.Here, we have considered the contribution to the far-field spatial decay of the net displacement in Eq. (81) from Eqs. (67)–(69) for the flow at order α2, Eq. (77) for the displacement of a tracer at order α, and Eqs. (37), (41) and (42) for the flow at order α. The result for the far-field net displacement at order α2 in Eq. (84) originates from the far-field source dipole in the translation contribution u2,0(T) to flow during a scan at order α2 (Eq. (70)), explained physically in Sect. 3.9.2.
Combining Eqs. (83) and (84), we therefore obtain the dimensional far-field net displacement of a tracer asΔX(t0)∼π-120αβ+112α2ΔT02U∫-t0t0A(t)2dta3R3(2cosθeR+sinθeθ)t=-t0.Hence, the far-field average velocity, over one scan, of the material point with initial position X0 is given byΔX(t0)2t0∼π-140αβ+124α2ΔT02Ut0∫-t0t0A(t)2dta3R3(2cosθeR+sinθeθ)t=-t0.This is a hydrodynamic source dipole in three dimensions. Its strength scales quadratically with the peak temperature change ΔT0 and linearly with the speed of heat-spot translation U (and hence linearly with frequency of scanning), while the direction depends on the thermal expansion coefficient α and thermal shear viscosity coefficient β purely through the factor of -140αβ+124α2. These two terms reflect and quantify the two different physical effects driving net transport: the interplay between thermal expansion and thermal shear viscosity changes at order αβ, and thermal expansion by itself at order α2. If the thermal expansion coefficient α and thermal shear viscosity coefficient β are both positive, as is the case for many liquids, then the two effects compete; the larger prefactor for the order-α2 term could potentially compensate for a value of α smaller than β. However, for a negative thermal expansion coefficient α, the two terms in the prefactor are instead of the same sign, thus reinforcing each other.
Discussion
In this article, we considered thermoviscous and thermal expansion-driven fluid flow in three-dimensional, unbounded fluid: a new geometry compared with previous work on viscous fluid confined between parallel plates [30]. We first examined heat transport. By solving numerically for the temperature field induced by a scanning heat source in the limit relevant to experiments [25, 27], we motivated a simplified model of the temperature perturbation to act as an input to our flow model. We then derived analytically the fluid flow and net transport due to the scanning heat spot, to quadratic order in the thermal expansion coefficient α and thermal shear viscosity coefficient β. Our model included bulk viscosity, a key new physical ingredient present because of the compressible nature of the flow. Bulk viscosity did not play a part in the parallel-plate setup due to the geometry, as shown using a scaling argument [30, 31]. Here, for the three-dimensional, unconfined fluid, we treated the bulk viscosity as a general function of space and time, without specifying its temperature dependence. Through our analysis, we found that even though the bulk viscosity impacts the pressure field, it does not influence the fluid velocity. We obtained the same physics and qualitative results as the parallel-plate setup. Specifically, the leading-order instantaneous flow during one scan is driven by thermal expansion via the spatio-temporally varying temperature field, occurring at order α and linear in the peak temperature change. In contrast with this, the leading-order net transport occurs at both order αβ and order α2, quadratic in the temperature change. The two typically competing contributions at these two orders originate from two fundamentally different physical mechanisms, with the order-αβ thermoviscous term generated by the interplay between thermal expansion and thermal shear viscosity changes, while the order-α2 effect arises purely from thermal expansion. The far-field average velocity of tracers is given by a hydrodynamic source dipole, but now in three dimensions.
We now comment on the validity of our spherical approximation for the heat spot. In existing microfluidic experiments, the scan path can have length on the order of 50μm, but there can be a relatively large length scale of heat absorption by the fluid, perpendicular to the scan path, of around 300μm. However, for a scan path long enough that heat spot is small in comparison, the geometry would be approximately axisymmetric; our model could therefore serve as a first approximation in this situation. Furthermore, methods of achieving thermoviscous transport in three dimensions that have been suggested include highly focused heating of the fluid [24].
Our work provides a first model of fully three-dimensional thermoviscous and thermal expansion-driven net flows, which we anticipate will be key to explaining experimental data (Moritz Kreysing, personal communication). The unconfined geometry of our model may help with understanding experiments where boundaries are far from the heat spot [25] or where the lubrication approximation (employed in theory for the parallel-plate setup [21, 30]) is no longer valid. Our theory will also allow us to analyse the effect of boundaries and confinement on three-dimensional thermoviscous and thermal expansion flows and transport in future modelling work. The quantitative theoretical results, in combination with scan-path selection via feedback algorithms [21, 22–23] or global optimisation [41], could contribute to the design of new experiments that apply net thermoviscous and thermal expansion flows to trap or manipulate particles in three dimensions.
Acknowledgements
We thank Moritz Kreysing for helpful discussions and feedback. We gratefully acknowledge funding from the Engineering and Physical Sciences Research Council (studentship to W.L.) and Trinity College, Cambridge (Rouse Ball and Eddington Research Funds travel grant to W.L.).
DeclarationsConflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.
ReferencesChengKGuoJFuYGuoJActive microparticle manipulation: recent advancesSens Actuators A202132210.1016/j.sna.2021.112616112616NilssonJEvanderMHammarströmBLaurellTReview of cell and particle trapping in microfluidic systemsAnal Chim Acta2009649214115710.1016/j.aca.2009.07.017ZhangZWangXLiuJDaiCSunYRobotic micromanipulation: fundamentals and applicationsAnnu Rev Control Robot Auton Syst20192118120310.1146/annurev-control-053018-023755MoffittJRChemlaYRSmithSBBustamanteCRecent advances in optical tweezersAnnu Rev Biochem200877120522810.1146/annurev.biochem.77.043007.090225Favre-BulleIAStilgoeABScottEKRubinsztein-DunlopHOptical trapping in vivo: theory, practice, and applicationsNanophoton2019861023104010.1515/nanoph-2019-0055PalimaDGlückstadJGearing up for optical microrobotics: micromanipulation and actuation of synthetic microstructures by optical forcesLaser Photon Rev20137447849410.1002/lpor.201200030De VlaminckIDekkerCRecent advances in magnetic tweezersAnnu Rev Biophys201241145347210.1146/annurev-biophys-122311-100544Baudoin M, Thomas J-L (2020) Acoustic tweezers for particle and fluid micromanipulation. Annu Rev Fluid Mech 52(1):205–234. https://doi.org/10.1146/annurev-fluid-010719-060154LiuSLinLSunH-BOpto-thermophoretic manipulationACS Nano20211545925594310.1021/acsnano.0c10427Lutz BR, Chen J, Schwartz DT (2006) Hydrodynamic tweezers: 1. Noncontact trapping of single cells using steady streaming microeddies. Anal Chem 78(15):5429–5435. https://doi.org/10.1021/ac060555yQianYNealeSLMarshJHMicroparticle manipulation using laser-induced thermophoresis and thermal convection flowSci Rep20201011916910.1038/s41598-020-76209-9WintererFMaierCMPernpeintnerCLohmüllerTOptofluidic transport and manipulation of plasmonic nanoparticles by thermocapillary convectionSoft Matter201814462863410.1039/C7SM01863KVaranakkottuSNGeorgeSDBaierTHardtSEwaldMBiesalskiMParticle manipulation based on optically controlled free surface hydrodynamicsAngew Chem Int Ed201352287291729510.1002/anie.201302111HuWFanQOhtaATAn opto-thermocapillary cell micromanipulatorLab Chip201313122285229110.1039/C3LC50389EArmaniMDChaudharySVProbstRShapiroBUsing feedback control of microflows to independently steer multiple particlesJ Microelectromech Syst200615494595610.1109/JMEMS.2006.878863TanyeriMJohnson-ChavarriaEMSchroederCMHydrodynamic trap for single particles and cellsAppl Phys Lett2010962222410110.1063/1.3431664ShenoyARaoCVSchroederCMStokes trap for multiplexed particle manipulation and assembly using fluidicsProc Natl Acad Sci USA2016113153976398110.1073/pnas.1525162113SchneiderTMMandreSBrennerMPAlgorithm for a microfluidic assembly linePhys Rev Lett2011106910.1103/PhysRevLett.106.094503094503TuMQNguyenHVFoleyEJacobsMISchroederCM3D manipulation and dynamics of soft materials in 3D flowsJ Rheol202367487787710.1122/8.0000600GonzalezJGopinathanALiuBSymmetry-group-protected microfluidicsPhys Rev Res20246210.1103/PhysRevResearch.6.023234023234ErbenESeelbinderBStoevIDKlykovSMaghelliNKreysingMFeedback-based positioning and diffusion suppression of particles via optical control of thermoviscous flowsOpt Express20212919302723028310.1364/OE.432935Stoev ID, Seelbinder B, Erben E, Maghelli N, Kreysing M (2021) Highly sensitive force measurements in an optically generated, harmonic hydrodynamic trap. eLight 1, 1–9 https://doi.org/10.1186/s43593-021-00007-7ErbenELiaoWMinopoliAMaghelliNLaugaEKreysingMOpto-fluidically multiplexed assembly and micro-roboticsLight Sci Appl20241315910.1038/s41377-024-01406-4WeinertFMBraunDOptically driven fluid flow along arbitrary microscale patterns using thermoviscous expansionJ Appl Phys20081041010.1063/1.3026526104701MittaschMGrossPNestlerMFritschAWIsermanCKarMMunderMVoigtAAlbertiSGrillSWKreysingMNon-invasive perturbations of intracellular flow reveal physical principles of cell organizationNat Cell Biol201820334435110.1038/s41556-017-0032-9SeelbinderBWagnerSJainMErbenEKlykovSStoevIDKrishnaswamyVRKreysingMProbe-free optical chromatin deformation and measurement of differential mechanical properties in the nucleuseLife2024137642110.7554/eLife.76421Minopoli A, Wagner S, Erben E, Liao W, Stoev ID, Lauga E, Kreysing M (2023) ISO-FLUCS: symmetrization of optofluidic manipulations in quasi-isothermal micro-environments. eLight 3(1):16. https://doi.org/10.1186/s43593-023-00049-zLuWGelfandVIGo with the flow - bulk transport by molecular motorsJ Cell Sci2023136510.1242/jcs.260300260300GoldsteinREMeentJ-WA physical perspective on cytoplasmic streamingInterface Focus2015542015003010.1098/rsfs.2015.0030LiaoWErbenEKreysingMLaugaETheoretical model of confined thermoviscous flows for artificial cytoplasmic streamingPhys Rev Fluids20238310.1103/PhysRevFluids.8.034202034202WeinertFMKrausJAFranoschTBraunDMicroscale fluid flow induced by thermoviscous expansion along a traveling wavePhys Rev Lett20081001610.1103/PhysRevLett.100.164501164501GosseCCroquetteVMagnetic tweezers: micromanipulation and force measurement at the molecular levelBiophys J 20028263314332910.1016/S0006-3495(02)75672-5HappelJBrennerHLow Reynolds number hydrodynamics1965Englewood CliffsPrentice HallHolmesMParkerNPoveyMTemperature dependence of bulk viscosity in water using acoustic spectroscopyJ Phys201126901201110.1088/1742-6596/269/1/012011TemkinSElements of acoustics1981New YorkWileyDukhin AS, Goetz PJ (2009) Bulk viscosity and compressibility measurement using acoustic spectroscopy. J Chem Phys 130(12):124519. https://doi.org/10.1063/1.3095471SlieWMDonforARJrLitovitzTAUltrasonic shear and longitudinal measurements in aqueous glycerolJ Chem Phys196644103712371810.1063/1.1726524XuJRenXGongWDaiRLiuDMeasurement of the bulk viscosity of liquid by Brillouin scatteringAppl Opt200342336704670910.1364/AO.42.006704DinçerİZamfirescuCDrying phenomena: theory and applications2016West SussexWileyYarivEBrennerHFlow animation by unsteady temperature fieldsPhys Fluids20041611959810.1063/1.1801091Erben E, Saraev I, Liao W, Nan F, Lauga E, Kreysing M (2025) Model-based thermoviscous micromanipulations via global optimization (submitted)Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.