2.1 Principles of the method

Figure 1 presents the principle of operation of the thermocapillary fluidic shaping method and the experimental set-up used for its implementation. When localized heating is imposed on a gas–liquid interface, the liquid will flow outward along the free surface from warmer regions with lower surface tension to colder regions with higher surface tension. This results in local deformation of the liquid film (Reference Tan, Bankoff and DavisTan, Bankoff, & Davis, 1990), and is known as the thermocapillary effect. The magnitude of the surface deformation is proportional to the local temperature gradients. We found that a convenient way to control these temperature gradients is to use light projection, whose spatial intensity can be easily programmed, to heat an array of absorbing metal pads patterned on the bottom of the fluidic chamber. The system uses a standard <3 W microscope light emitting diode (LED) source and a standard DMD. The resulting increase in temperature is of the order of few degrees or less, and is sufficient to drive significant deformations.

The spatial resolution of the thermal patterning depends on the material properties of the substrate, the resolution of the projected pattern and the geometry of the light-absorbing metallic pad array. The composition of the metal pads was designed to maximize light absorbance across the ultraviolet-visible (UV-Vis) spectrum of the light source, as shown in figure 1(c). Section 2 in the supplementary movie and material available at https://doi.org/10.1017/flo.2022.17 provides characterization of the deformation magnitude and resolution in our system as a function of the illumination intensity and pad array density. As expected, the magnitude of the resulting deformation is proportional to both the illumination intensity and the area fraction occupied by the pads. The precise deformation magnitude and resolution depend on the specific heating pattern, with some tradeoff between the two (see the theory section). However, typical values are of the order of 3 μm maximum deformation and 5 μm mm⁻¹ spatial resolution.

2.2 The inverse problem − theory and implementation

Consider an open microfluidic chamber of length $\mathrm{{\cal L}}$, containing a thin liquid film of thickness ${h_0}$, as depicted in figure 1(a). The chamber is heated from below by a light pattern that is absorbed by the metal pads array, resulting in temperature variations over a characteristic length l, ${h_0} \ll l \ll \mathrm{{\cal L}}$. For small variations in temperature, we can assume that the surface tension, $\sigma $, changes linearly with temperature, $\theta $, according to

(2.1)\begin{equation}\sigma = {\sigma _r} - \beta (\theta - {\theta _r}),\end{equation}

where $\beta \equiv \partial \sigma /\partial \theta $ is defined as the thermocapillary coefficient and is assumed to be constant, ${\theta _r}$ is a reference temperature and ${\sigma _r}$ is the surface tension at that temperature. We assume that variations in temperature are sufficiently small such that changes in density and viscosity are negligible.

In supplementary material § 3 we provide a complete derivation of the governing equations. Briefly, we consider the continuity equation, the steady-state incompressible Navier–Stokes equations and the steady-state advection–conduction heat equation. At the bottom of the film $(z = 0)$ we require no-slip and no-penetration conditions on the flow, and set a fixed temperature distribution, ${\theta _b}$. At the liquid--air interface $(z = h(x,y))$ we impose the kinematic condition relating flow velocities to the shape of the deformed interface, a stress balance relating the stress tensor in the liquid to surface tension and heat loss through Newton's cooling law. Using these equations and boundary conditions, we derive the dimensionless steady-state equation of the deformed liquid film driven by a non-uniform temperature distribution

(2.2)\begin{equation}{\boldsymbol{\nabla} _2} \boldsymbol{\cdot} \left[ {\frac{S}{3}{H^3}{\boldsymbol{\nabla}_2}\nabla_2^2H - \frac{G}{3}{H^3}{\boldsymbol{\nabla}_2}H - \frac{1}{2}{H^2}{\boldsymbol{\nabla}_2}\left( {\frac{{{\vartheta_b} + \varTheta }}{{1 + BH}}} \right)} \right] = 0.\end{equation}

Here, $[X,Y,Z,H] = [x/l,y/l,z/{h_0},h/{h_0}]$ are the dimensionless coordinates and film thickness scaled by the characteristic lateral and vertical dimensions, and ${\boldsymbol{\nabla} _2} = (\partial /\partial X,\partial /\partial Y)$ and $\nabla _2^2 = ({\partial ^2}/\partial {X^2} + {\partial ^2}/\partial {Y^2})$ are the two-dimensional gradient and Laplacian operators. We scale the velocity as $[U,V,W] = (1/{U_0})[u,v,(1/\varepsilon )w]$, where $U_0 = \varepsilon \mathrm{\Delta }\sigma/\mu $ is the characteristic velocity as derived from the tangential stress balance, with $\Delta \sigma = \beta \Delta \theta $ being the characteristic variation in surface tension, and $\mathrm{\Delta }\theta = {\theta _{{b_{max}}}} - {\theta _{{b_{min}}}}$ the difference between the maximum and minimum temperatures at that surface. Also, $\mu $ is the dynamic viscosity of the liquid, and $\varepsilon = {h_0}/l \ll 1$ is the ratio of film thickness to characteristic lateral scale. We scale the temperature as $\vartheta = (\theta - {\bar{\theta }_b})/\mathrm{\Delta }\theta $, where ${\bar{\theta }_b}\; $is the average temperature of the bottom surface, $S = {\varepsilon ^3}({\sigma _r}/\mu {U_0}) = {\varepsilon ^2}({\sigma _r}/\mathrm{\Delta }\sigma )$ is the surface tension number $B = {{\mathscr{h}}h_{0}}/{k}$ is the Biot number, ${\mathscr{h}}$ is the heat transfer coefficient at the liquid–air interface, k is the thermal conductivity, $\varTheta = ({\bar{\theta }_b} - {\theta _\infty })/\mathrm{\Delta }\theta $ is a ratio of temperature differences, where ${\theta _\infty }$ is the ambient temperature, and $G = \varepsilon (\rho h_0^2/\mu {U_0})g = (\rho h_0^2/\mathrm{\Delta }\sigma )g$ is the dimensionless gravity.

Subjected to appropriate boundary conditions on the thickness of the liquid film (e.g. far from the heating region, the film thickness can be assumed to be known and uniform), (2.2) can be solved numerically to obtain a solution to the direct problem − the spatial deformation resulting from any prescribed temperature distribution ${\vartheta _b}(x,y)$. However, for engineering purposes, the inverse problem − seeking the temperature map that would yield a desired deformation − is of greater value.

We approach the inverse problem by first considering the one-dimensional case described by

(2.3)\begin{equation}\; \; \; {\partial _\textrm{X}}\left[ {\frac{S}{3}{H^3}\partial_X^3H - \frac{G}{3}{H^3}{\partial_\textrm{X}}H - \frac{1}{2}{H^2}{\partial_X}\left( {\frac{{{\vartheta_b} + \varTheta }}{{1 + BH}}} \right)} \right] = 0,\end{equation}

where X is the non-dimensional spatial coordinate. We assume that the temperature variations occur far from the boundaries, $X ={\pm} \mathrm{{\cal L}}$, and therefore at the boundaries the film can be assumed to be flat

(2.4a–d)\begin{equation}H({\pm} \mathrm{{\cal L}}) = 1,\quad \partial _X^nH({\pm} \mathrm{{\cal L}}) = 0\{ n = 1,2,3\} ,\quad {\vartheta _b}({\pm} \mathrm{{\cal L}}) = {-} \varTheta ,\quad {\partial _X}{\vartheta _b}({\pm} \mathrm{{\cal L}}) = 0.\end{equation}

Integrating equation (2.3) in X and dividing by ${H^2}$ we obtain

(2.5)\begin{equation}\frac{S}{3}H\partial _X^3H - \frac{G}{3}H{\partial _X}H - \frac{1}{2}{\partial _X}\left( {\frac{{{\vartheta_b} + \varTheta }}{{1 + BH}}} \right) = \frac{{{C_0}}}{{{H^2}}},\end{equation}

where ${C_0}$ is an integration constant. The left-hand side can be integrated by parts, and vanishes when using the boundary conditions (2.4a–d)

(2.6)\begin{equation}\int_{ - \mathrm{{\cal L}}}^\mathrm{{\cal L}} {\left[ {\frac{S}{3}H\partial_X^3H - \frac{G}{3}H{\partial_\textrm{X}}H - \frac{1}{2}{\partial_X}\left( {\frac{{{\vartheta_b} + \varTheta }}{{1 + BH}}} \right)} \right]\,\textrm{d}X = 0} .\; \end{equation}

Since the integration of the right-hand side of (S4) is non-zero, the coefficient ${C_0}$ must vanish, ${C_0} = 0$.

The one-dimensional steady state (2.5) can thus be simplified to

(2.7)\begin{equation}\frac{S}{3}H\partial _X^3H - \frac{G}{3}H{\partial _X}H - \frac{1}{2}{\partial _X}\left( {\frac{{{\vartheta_b} + \varTheta }}{{1 + BH}}} \right) = 0.\end{equation}

Integrating equation (2.7) to an arbitrary position $\xi $ yields an explicit expression for the temperature variation, ${\theta _b}(\xi )$, required in order to create the desired deformation $H(\xi )$:

(2.8)\begin{equation}{\vartheta _b}(\xi ) = \frac{{1 + BH(\xi )}}{3}[2SH(\xi )\partial _\xi ^2H(\xi ) - S{({\partial _\xi }H(\xi ))^2} - G(H{(\xi )^2} - 1)] - \varTheta .\end{equation}

The two-dimensional case cannot be directly integrated. We therefore propose an ansatz based on the one-dimensional solution

(2.9)\begin{equation}{\vartheta _b} = \frac{{1 + BH}}{3}[2SH{\nabla ^2}H - S{(\boldsymbol{\nabla }H)^2} - G({H^2} - 1)] - \varTheta .\end{equation}

As detailed in the supplementary movie and material, by substituting (2.9) into the governing equation (2.2), we obtain the difference between the approximate (ansatz) solution and the exact solution

(2.10)\begin{equation}\boldsymbol{\epsilon } = \frac{S}{3}[({\partial _Y}H\partial _{XY}^2H - {\partial _X}H\partial _Y^2H),({\partial _X}H\partial _{XY}^2H - {\partial _Y}H\partial _X^2H)].\end{equation}

The ansatz solution is thus exact when the non-normalized Gaussian curvature of the surface, $K = \partial _X^2H\partial _Y^2H - {(\partial _{XY}^2H)^2}$, vanishes. For other surfaces, the error would depend on the magnitude of the Gaussian curvature, and we provide a numerical assessment for this error in supplementary material § 3.

Figure 2 presents the use of the inverse problem for fabricating a desired topography. Figure 2(a) presents the prescribed topography composed of a set of concentric rings, with a maximum deformation relative to the baseline film thickness of $d_{desired} = 1.5\;\mathrm{\mu }\textrm{m}$. This topography is used as input into (2.3), with the constant parameters as defined in supplementary material § 4, resulting in a temperature distribution ${\vartheta _b}(x,y)$. Since we do not have a direct relation between the desired temperature and the required illumination intensity, we assume a linear dependence and perform a calibration measurement in which the illumination intensity map is set to $I(x,y) = {I_0}\; {\vartheta _b}(x,y)/\; \textrm{max}\{ {\vartheta _b}(x,y)\} $, where ${I_0}$ is the maximum intensity of the light source. The resulting intensity map is presented in figure 2(b), and an image of the actual projected intensity map is presented in figure 2(c). We measure the resulting maximum deformation, ${d_{max}}$ (which for this pattern was 2.2 μm), and repeat the projection using an intensity map that is scaled by ${d_{desired}}/{d_{max}}$ (i.e. 1.5/2.2). Figure 2(d) presents the resulting topography of the (solidified) deformed liquid film. There is a good qualitative agreement between the desired pattern and the resulting one yet deviations are also clearly visible. We further discuss these deviations in the Conclusions section.

Figure 2. Demonstration of a complete workflow for an inverse-problem solution. (a) The desired topography defined in terms of the variation from a reference height. (b) The intensity map required in order to produce the desired deformation, as obtained from solving the inverse problem. (c) Image of the illumination pattern as obtained on the bottom surface of the fluidic device. (d) The resulting topography after curing the polymer as measured by the DHM (Reference Cuche, Marquet and DepeursingeCuche, Marquet, & Depeursinge, 1999).

2.3 Applications

Figures 3–5 present several practical diffractive elements produced using the thermocapillary fluidic shaping approach. Diffraction gratings are optical elements with a periodic structure, that diffract the incident light into several beams traveling in different directions, depending on the wavelength and the spatial periodic structure defining the angular spread. They are commonly used in monochromators for spectroscopic instruments (Reference Loewen and PopovLoewen & Popov, 2018).

Figure 3 presents a linear diffraction grating with a peak-to-peak periodicity of 400 μm produced by projection of a sinusoidal illumination pattern. Figure 3(a) presents the projected pattern side by side with the resulting deformation, and figure 3(b) presents the magnitude of the deformation along the centreline. We test the resulting element by measuring the diffractive angle for several wavelengths produced by monochromatic LEDs. Since the LEDs have a spectral bandwidth of approximately 50 nm, we use a weighted average of their spectrum (as provided by the manufacturer) to define a representative wavelength for the analysis. As shown in figure 3(c), the results are in excellent agreement with the expected angle set by the element's periodicity. Figure 3(d) also presents a visual demonstration of white light diffracted through the grating and directly imaged by a camera sensor, clearly showing the separation of bands.

Figure 3. Fabrication and testing of a linear diffraction grating. (a) Side-by-side image showing the projected light pattern (left) and the measured resulting topography (right). Regions with a higher projected intensity, and therefore a higher temperature, correspond to valleys in the topography. (b) The deformation magnitude along the centre of the element, showing the spatial frequency (400 $\mu$m) and the resulting amplitude (0.4 $\mu$m) of the grating. (c) The measured diffractive angle for three wavelengths passing through the grating. The reciprocal of the slope of a least-square line matches well the designed grating spatial frequency. (d) Image of white light diffraction. The numbers along the scale indicate the centre of each band along a horizonal centreline.

The depth of field of any optical system can be extended by reducing its numerical aperture, yet at the cost of reducing the amount of light collected. A number of studies have shown the ability to overcome this limitation and extend the depth of field without loss of photons, by positioning an appropriate phase mask in the systems’ aperture (Reference Ben-Eliezer, Marom, Konforti and ZalevskyBen-Eliezer, Marom, Konforti, & Zalevsky, 2005; Reference Dowski and CatheyDowski & Cathey, 1995). Figure 4 presents the fabrication and testing of an extended depth of focus (EDOF) phase mask designed by formulating a phase retrieval task as described by Nehme et al. (Reference Nehme, Ferdman, Weiss, Naor, Freedman, Michaeli and ShechtmanNehme et al., 2021). The desired mask mould is translated into a light pattern (figure 4b) using the inverse-problem solution, and the thermocapillary fluidic shaping method is used to create a solid mould. We then cast Polydimethylsiloxane (PDMS) onto the mould to create the phase mask itself, shown in figure 4(c). We position the mask in the Fourier plane of a 4f system and use it to image three targets positioned at different distances from the camera. As illustrated in figure 4(d), to ensure that the presence of the mask itself does not reduce the numerical aperture, we use a smaller fixed aperture in close proximity to the mask, acting as the aperture stop of the system. Figures 4(e) and 4(f) clearly show the additional depth of field gained by using the mask, yet with some loss of contrast.

Figure 4. Fabrication and demonstration of an extended depth of field phase mask. (a) The design of the EDOF mask based on Nehme et al. (Reference Nehme, Ferdman, Weiss, Naor, Freedman, Michaeli and ShechtmanNehme et al., 2021). (b) The projection pattern for a negative mould, as solved by the inverse problem. (c) The obtained phase mask as measured by the DHM. (d) Schematic illustration of the optical set-up. We position the phase mask in the Fourier plane of a 4f system, and image three targets positioned 2.5 cm from one another. We used an additional small fixed aperture to ensure that the phase mask itself does not reduce the aperture stop. (e) The resulting image without the phase mask, showing only the centre target in focus. (f) The resulting image with the EDOF phase mask showing all three targets in focus.

Three-dimensional localization microscopy uses phase masks positioned in the Fourier plane of a microscope to modify the PSF in order to encode three-dimensional information (Reference Pavani, Thompson, Biteen, Lord, Liu, Twieg and MoernerPavani et al., 2009). Figure 5 presents the use of the inverse-problem approach for fabrication of such a phase mask. Figure 5(a) presents the topography of the desired ‘tetrapod’ mask (Reference Shechtman, Sahl, Backer and MoernerShechtman et al., 2014), designed to provide depth information over a range of ±2 μm. Here, too, we use the desired mould topography as input to (2.3) and obtain the intensity map presented in figure 5(b). The parameters used in the solution are identical to those that were used in figure 4 and that are detailed in supplementary material § 3. Figure 5(c) presents the topography of the resulting phase mask created by casing PDMS onto the resulting mould, as measured by the digital holographic microscope (DHM). The first column in figure 5(d) presents the measured PSF of a single emitter (a 200 nm fluorescent bead) at different locations along the optical axis, in the absence of any phase mask. The second column presents the measured PSF of the same emitter, when using a phase mask produced using a 3-step lithography process. The last column in figure 5(d) presents the measured PSF as obtained using our phase mask fabricated by the thermocapillary fluidic shaping method. Figure 5(e) demonstrates the use of this mask for three-dimensional tracking of three 200 nm beads undergoing Brownian motion (see supplementary movie 1 and material) in a liquid environment, and of one stationary bead (attached to the bottom surface).

Figure 5. Fabrication and demonstration of a saddle-point phase mask for three-dimensional localization microscopy. (a) Topography map of a saddle-point phase mask microfabricated using standard lithography and ion etching processes. (b) The illumination pattern required for deforming the liquid film into a negative mould for the phase mask, as obtained from the inverse-problem solution. (c) The DHM measurement of the resulting solidified PDMS mask, cast on the mould fabricated by thermocapillary fluidic shaping. (d) The PSF of a 200 nm diameter fluorescent microsphere at various z positions, as obtained from (i) a measured standard PSF with no phase mask (ii) a measured PSF using a lithography phase mask, (iii) a measured PSF using the mask obtained through thermocapillary fluidic shaping. (e) Tracking of three nanobeads undergoing Brownian motion and one stationary bead (in green) in a liquid, using the fabricated mask.