Experiments

The first appearance of the microscopic cavitation domain was detected at Er≃200. As shown in Fig. 1a,b, the vapour phase localizes at the hydrodynamic stagnation point downstream of the micron-sized pillar. As the streamlines of the flowing 5CB divide upon approaching the micropillar, both Er and Re first increase locally, reach maximum values at the point of minimum separation and reattain their original values as the flows converge downstream of the pillar. With an average flow speed of 800 μm s−1 (Er=200), the velocity reaches a maximum of 4,000 μm s−1 at the constriction (minimum distance of 10 μm between the pillar surface and the channel wall). This corresponds to Er=945 and Re=0.1, a hydrodynamic Stokes regime. The resulting drop in the hydrodynamic pressure corresponds to ≃8 mPa, corroborating that the absolute pressure decreases below the vapour pressure of nematic 5CB (≃0.1 mPa at 298 K47).

Figure 1: Cavitation in nematic LC due to flow past a micron-sized pillar. (a) Polarized optical micrograph and (b) white light micrograph show steady-state cavitation domain in nematic LC, 5CB. The channel and pillar surfaces have homeotropic boundary condition. Right panel: Magnified projection of the cavitating domain. (c) Minimum intensity projection of the scattered light over time due to flowing disclination lines. The disclinations were used to trace the streamlines of the LC flow past the micropillar. The cavitating domain is locally stabilized at the hydrodynamic stagnation point (indicated by the red arrowhead) downstream of the micron-sized pillar. For the micrographs shown here, the pillar diameter is 2r=80 μm, placed within a 100 μm wide and 15 μm deep microchannel. (d) Map of the local nematic order parameter in proximity of the cavitation domain obtained from MD simulations (see attached colour bar). The area of solid red colour corresponding to vanishing nematic order represents the cavitation. Small dashes in the maps represent the local director field. The length of a dash represents its three-dimensional orientation: shorter dashes have an orientation closer to the normal to the plane of the map. (e) Magnified view of the cavitating domain observed between crossed polarizers. The intensity of the transmitted light (normalized by the maximum intensity) is measured along the dashed lines: red, blue and green, and plotted in (f). For each line, the maximum value of the normalized intensity (bright regions) is recorded outside the cavitation domain (bulk director is oriented at 45° relative to the polarizers). At the apex of the cavitation domain and further downstream, the director is oriented parallel to the flow direction that appears dark between the crossed polarizers. The gas-filled cavitation domain also appears dark (minimum intensity) due to the total extinction of transmitted light. (g) Cavitation domain was also observed under planar surface anchoring conditions, visualized here between crossed polarizers. Full size image

Figure 1c shows the downstream hydrodynamic stagnation point (red arrowhead), visualized using a differential intensity projection of the light scattered by the flowing disclinations (Methods). Due to the high Er flows, singular topological defects readily fills up the nematic bulk that scatter white light due to the lensing effect48. These freely flowing disclinations (both ±1/2 and ±1 defects are present) follow the streamlines, and accurately map the laminar flow field around the micropillar. Between the crossed polarizers, the cavitation domain has a dark appearance due to complete extinction of the transmitted white light. As shown in Fig. 1b, the interface between the vapour and liquid phases of 5CB is distinctly visible under white light. Cavitation was also observed within microfluidic devices possessing planar surface anchoring (Fig. 1g). Irrespective of the surface anchoring, no cavitation was observed in nematic flows at Er<200. In all cases, the nematic director field, average orientation of the LC molecules, visualized using polarization microscopy, was observed to be aligned by the flow field (nematic 5CB is a flow aligning LC32). This agrees very well with flow-induced director field described previously49.

Figure 1d maps the local nematic order in proximity of the cavitation domain, obtained through MD simulations. In agreement with the polarized optical measurements (Fig. 1e), our simulations predict that the LC molecules anchor homeotropically at the fluid–vapour interface. As shown in Fig. 1f, the normalized intensity of the transmitted light, plotted for three adjacent sections on the cavitation domain (red, blue and green dashed lines), attains maximum value just outside the cavitation domain. The maximum intensity corresponds to bulk director orientation of 45° relative to the crossed polarizers. Close to the apex of the cavitation domain, the director is oriented along the flow direction. This appears dark between crossed polarizers as the director is parallel to one of the polarizers. A minimum intensity is also observed inside the cavitation domain due to total extinction of the white light passing through the gas phase between crossed polarizers. The intermediate values of the intensity at the cavitation interface suggest that the molecules have a weak homeotropic alignment at the interface, as also confirmed by our MD simulations. Measurements on 5CB show that the molecules align homeotropically at the fluid–vapour interface50,51,52, in agreement with existing studies of the LC–air interface for 5CB53. Our MD simulations correctly capture this behaviour that lends confidence to the predictions about the underlying physics. The disclination lines found in the experiments are not present in the MD system (Fig. 1d) as they require much larger system sizes to develop.

Upon appearance, the cavitation volume increases over time, and typically within few hours (4–7 h) of steady flow of the nematic 5CB, saturates into a constant size. Figure 2a,b show polarization optical microscopy (POM) micrographs of a cavitating domain developing over 5 h at . The time taken to attain the steady-state volume varies in experiments and, in addition to the roughness of the pillar surface (higher roughness promotes cavitation), it showed a qualitative dependence on the depth of the microchannel. At a given Er, it took longer to reach the saturated volume in a shallow channel (d<12 μm), compared with that in deeper channels. However, once formed, the cavitation domain in a shallow channel was observed to be stable for longer times, possibly due to strong pinning of the vapour–liquid interface at the top and bottom surfaces of the microchannel. The POM micrographs with λ-plate (Fig. 2b) confirm the absence of the nematic phase within the cavitation domain. While cavitation can be initiated by pre-existing gas bubbles entrapped in the asperities of the confining solid surfaces, we do not expect this to play a determining role in our experiments. Under similar surface conditions, no cavitation was observed when 5CB in the isotropic phase was flowed, and in experiments with deionized water. The inception of cavitation only in the nematic 5CB suggests that cavitation is triggered due to the complex anisotropic nature of the flowing material. Cavitation was also reproduced with obstacles having noncircular cross-sections, including semicircular and square geometries (Fig. 2c). At high flow speeds (Er>500), the domain volume was observed to reduce in size. This is counterintuitive, since high flow speed (high Re) is known to favour the onset of cavitation in isotropic fluids. However, in the present setting, we speculate that the strong viscous forces disrupt the pinned vapour–liquid interface, and shear off minuscule gas bubbles from the saturated cavitation domain. This can reduce the domain volume, triggering a dynamic tradeoff between the growth and decay processes that regulate the final cavitation volume at a given Er.

Figure 2: Cavitation domain growth. (a) Polarized optical micrographs represent the growth of the cavitation volume over time. The total time elapsed is 5 h. The minimum intensity regions (dark appearance) are observed either due to the extinction of the transmitted light as it passes through the vapour phase (in the cavitation domain), or through the bulk nematic phase aligned parallel to one of the polarizers. The two cases are distinguished using a λ-plate. (b) Introduction of the λ-plate confirms the absence of the nematic phase in the cavitation domain, and distinguishes it from the bulk nematic aligned parallel to the polarizer. (c) Cavitating domains were observed upon LC flow past different obstacle geometries: semicircle (left) and square (middle); and at different channel depths (right, d≈10 μm). Shallow channels require high Er numbers for cavitation to be triggered. Full size image

Simulations

To obtain an in-depth understanding of the physical processes that underlie the cavitation phenomenon, we turn to a theoretical modelling of our experimental setup. Due to the limitations in computer performance, the length scales accessible to MD simulations are still too small compared with the micron-sized cavitation domain in our experiments. However, the MD simulations do provide us access to the early stages of the LC cavitation, in addition to the structural information that is not easily tractable with phenomenological theories of LCs (see Methods and Supplementary Information). To characterize the onset and evolution of cavitation within the LC fluid, in Fig. 3a–c we present the maps of the local density in a x–y cross-section of the system. For the system at rest, that is, without flow, the density is homogeneous throughout the system. For moderate flow conditions the density remains homogeneous across the system (see Fig. 3a). However, at high flow rates, starting at Er≃897, the density across the system starts to change. Figure 3c shows a well-developed cavitation domain at Er≃1,374. Specifically, we observe a drop in the density in the downstream wake of the pillar (Fig. 3b). The density of this area corresponds to that of vapour phase, and therefore demonstrates hydrodynamic cavitation. The size of the vapour domain first increases till Er≃1,374, and then grows at a reduced rate. At Er≃1,586 the cavitation domain reaches the maximum size and does not grow further.

Figure 3: MD simulation of the hydrodynamic states. (a–c) Maps of the magnitude of local density in the x–y plane located at z=0 (see attached color bar). The circle represents the cylindrical pillar and the arrow above each column indicates the direction of flow, at different Ericksen numbers: (a) Er≃65, (b) Er≃897 and (c) Er≃1,374. (d–f) MD simulated local pressure plots and (g–i) corresponding local velocity fields, carried out at the same Er as in (a–c). Colour bars attached to plots in the top panel refer to all plots in the same column. Full size image

We notice that the Er at which cavitation initiates is about a factor of four higher than in our experiments. This apparent discrepancy is resolved by estimating the nematodynamic parameters locally around the micron-sized obstacle, rather than the far-field values. Locally, the Er reaches much higher values due to the increased flow speed through the constriction between the pillar and channel walls44. When cavitation first occurs in our experiments, at a far-field Er of 200 (v=800 μm s−1), the corresponding local Ericksen number Er loc =945 (v loc =4,000 μm s−1). This value is intriguingly close to Er=897 obtained in the MD simulations. We note that the local Reynolds number still falls within the Stokes regime . The close agreement between the experiments and simulations prompts us to conclude that the microfluidic experiments and the MD simulations represent equivalent hydrodynamic conditions. Thus, despite the difference in length scales, by analysing the system over comparable nematodynamic conditions, we are able to capture the underlying mechanism for cavitation under relevant boundary conditions.

The development of cavitation depends on the local pressure landscape. Figure 3d–f illustrate maps of the local pressure in the x–y cross-section of the system. The system at rest exhibits a homogeneous local pressure, as one would expect(not shown here). From the map in Fig. 3d it is clear that this applies to small Ericksen numbers Er≃65 as well. However, as shown in Fig. 3e, for Er≃897 (the onset of cavitation), this situation changes drastically. On the upstream side, the local pressure builds up because of the flow obstruction by the pillar. This is accompanied by a pressure drop on the downstream side, just behind the cylindrical pillar. After passing the constriction the fluid expands freely, and thus resulting in the drop in local pressure. This effect becomes even more pronounced if flow is further increased as shown in Fig. 3f. Notice that within the cavitation area the local pressure equals to the pressure of an ideal gas.

Figure 4 presents the corresponding evolution of the local pressure in the vertical plane (y–z plane) at specific distances behind the cylindrical pillar. For reasons of symmetry, and to obtain improved statistics, the values of the local pressure have been averaged along the channel depth (z-axis), and are presented across the channel width (y-axis). To avoid confinement effects, we exclude the regions close to the planar substrates from our calculation. For small flow velocities, the local pressure behind the cylindrical pillar is, on average, constant along the y-axis. If flow is increased (see Fig. 4a), the local pressure starts decreasing behind the cylindrical pillar. At the threshold value for cavitation, Er≃786, the local pressure attains negative values. Hence, molecules in this area are pushed out. This occurs just before a cavitation volume nucleates. As soon as a cavitation volume is visible (Er≃897) the local pressure in that area is equal to the ideal gas pressure.

Figure 4: Mapping the local pressure. (a) Dependence of the local pressure at distance 0.7 behind the surface of the cylindrical pillar in the direction of flow, plotted across the channel width (y-axis) for Er≃65 (red line), Er≃583 (brown dashed line), Er≃786 (orange dashed line) and Er≃897 (cyan dotted line). The inset in (a) shows a sketch of different positions behind the cylindrical pillar at which the local pressure was measured. (b–e) Plots of local pressure at specific distances behind the surface of the cylindrical pillar for Er≃65: (b) at a distance of 0.7 behind the pillar, (c) at 1.3, (d) at 1.7 and (e) at 2.7. Full size image

The oscillations visible in Fig. 4a for Er≃65 are due to the layering effects54,55 of the molecules close to the pillar. This leads to changes in pressure in this region located at the centre of the y-axis. If the local pressure is calculated further away from the surface of the cylindrical pillar, layering effects become less prominent, and as shown in Fig. 4b–e, the local pressure minimum (at y=0) vanishes. From a molecular standpoint, the local pressure oscillations close to the surface of the cylindrical pillar support the nucleation of a cavitation domain. This is consistent with the fact that hydrodynamic cavitation requires a nucleation site such as a surface.

Figure 3g–i maps the local velocity over the mid-plane cross section (x-y plane) of the nematofluidic system. At small flow rates, Er≃65, the local velocity is rather low, as indicated by the random orientations of the arrows representing the flow direction (Fig. 3g). At the onset of cavitation, Er≃897, a net flow is clearly present (Fig. 3h). The local velocity in the constriction is considerably higher than that in front of and behind the cylindrical pillar. The relative velocity difference becomes stronger at still higher flows. Figure 3i presents the velocity map at Er≃1,374 in the vicinity of the obstacle and the cavitation domain. Due to lack of sufficient statistics within the gaseous phase, the local convective velocity was set to zero in the cavitation area for calculating the resulting velocity map.

We apply the density maps presented in Fig. 3a–c for estimating the volume, V c , of the cavitation domain. The total volume was obtained by adding up all the cubic volumes for which the local density ρ l ≤0.15. Figure 5 plots the variation of V c with the Reynolds number Re, where Re cr denotes the critical Reynolds number for the cavitation inception. Alongside anisotropic system described so far, we have also performed MD simulations for the same LC model in the isotropic phase (T=1.1). The two systems differ in their viscosity η: the nematic LC (η≃37) is more viscous than the isotropic LC (η≃21). Thus, from the definition of Re, it follows that the nematic LC cavitates at lower critical Reynolds number Re cr than the isotropic LC. That fluids with lower viscosity require a larger Re cr to cavitate is consistent with our experiments. Neither the deionized water nor 5CB in isotropic phase elicited cavitation at Reynolds numbers in which the nematic 5CB cavitated. It is important to note here that, in our simulations, the Re cr for the nematic LC is of the same order of magnitude as the local Re in the microfluidic experiments (≃0.1).

Figure 5: Growth of cavitation volume. Cavitation volume plotted against Reynolds number Re for the nematic LC at T=0.88 (black triangle), at T=0.90 (red circle) and at T=0.92 (brown diamond), and for an isotropic LC (blue square). The inset shows the critical Reynolds number, Re cr , plotted against the average nematic order parameter , spanning temperatures from T=0.92 down to T=0.87. Full size image

For both the nematic and the isotropic phases of liquid crystals, we observe a steep growth of the cavitation domain after its formation, followed by a saturation to a plateau. The height of the saturation plateau is due to the finite size of the simulation box and, therefore, is system size dependent. However, the plateau value of the cavitation domain depends on the nature of the fluid, and can be ascribed to the differences in the liquid matrix structure, for example, different ordering among the three fluids. Figure 5 shows the growth of cavitation domains for the nematic LC at T=0.88 and T=0.92 that correspond to higher (η≃40) and lower (η≃34) viscosities, respectively, than what we have considered till now. While the qualitative trend is still retained, quantitative differences emerge, specifically in the growth and saturation of the cavitation domain. The lower the temperature (higher viscosity), the larger was the resulting cavitation volume. In addition to viscosity, the nematic order parameter also varies with T. The inset of Fig. 5 shows the dependence of Re cr on the average nematic order parameter for temperatures ranging from T=0.92 down to T=0.87. For values of in the nematic phase the critical Reynolds number Re cr decreases linearly with , demonstrating that cavitation is enhanced by stronger nematic alignment. This qualitative measurement serves as a first step towards understanding the influence of long range ordering on cavitation at micro scales.

Figure 6 shows the cavitation volume V c in terms of the inverse Euler number , a dimensionless number that is independent of the viscosity of the system. The ratio between the pressure gradient and the kinetic energy per volume is decisive for the development of a cavitating volume, as captured distinctly in Fig. 6. Independent of the model system cavitation occurs at the same Euler number . When cavitation first occurs in the experiments we find that is within similar order of magnitude as the simulations.

Figure 6: Cavitation volume and Euler number. Dependence of the cavitation volume on the inverse Euler number for the nematic LC (red circle) and isotropic LC (blue square). Full size image

Figure 7 maps the local nematic order parameter over the range of flow regimes considered here. The flow field reorients the nematic director from homeotropic (on the surfaces) to flow-aligned orientation (in the flowing matrix). Topological defects arise in the director field, close to the top and the bottom walls, where the cylindrical pillar intersects the channel surfaces (Fig. 7a). The singular defect loops are formed due to the homeotropic anchoring, both on the channel surfaces and on the pillar, and are consistent with the defect topology discussed in ref. 49. It is evident from Fig. 7d that over the x–y plane, located at the channel half depth (z=0), no defect is visible. The director field remains stable for small Er. However, for Er≥675 a single loop around the pillar stabilizes in the x–y mid-plane. The loop is deformed and extended towards the downstream direction along with the flow, shown in Fig. 7b,e. Additionally, there is a growing region of flow alignment behind the cylindrical pillar in the downstream direction. Upon increasing Er one can see that the loop becomes stretched further downstream (see Fig. 7c). However, the overall defect topology, especially downstream behind the pillar, is increasingly smeared, possibly due to the appearance of the vapour phase at high Ericksen numbers (Fig. 7f). Changing the surface anchoring from homeotropic to planar did not produce any qualitative change in our results. This agrees well with the experiments, where too we have observed that cavitation in nematic phase was independent of the nature of the surface anchoring.

Figure 7: Flow-induced local nematic order. (a–c) Magnitude of the local nematic order in the x–z plane located at y=0 (see attached colour bar). The cylindrical grey area represents the cylindrical pillar and the arrow above each plot indicates the direction of flow. Small dashes in the maps represent the local director field such that the length of a dash denotes its three-dimensional orientation: shorter dashes have an orientation closer to the normal to the plane of the map. The Ericksen numbers considered are (a) Er≃130, (b) Er≃675 and (c) Er≃786. (d–f) Corresponding director fields over the x–y plane located at z=0. Full size image

Theoretical analysis

We theoretically study the cavitation in an anisotropic matrix and predict the role of nematic order parameter on cavitation inception. As the simplest case, let us consider a spherical environment containing the nematic liquid crystal (see Methods and Supplementary Fig. 1). The fluid is subject to a pressure drop ΔP. Under appropriate hydrodynamic conditions, a stationary pressure drop leads to cavitation in the nematic bulk. The first region of interest, the spherical cavitation domain with radius r c , encloses the vapour phase of the LC. From our foregoing experiments and computer simulations we know that the nematic LC molecules align homeotropically to the surface of the vapour phase. This we consider as the region 2: a spherical shell of thickness r d −r c with perfect homeotropic alignment to region 1. Finally, the rest of the system is the bulk nematic phase that we will treat in the mean-field approximation, namely, with a homogeneous nematic order parameter S. Locally, the order parameter may vary as a result of the inception and growth of the cavitation domain. The third region has a radius r L >r d , and excludes the sphere of radius r d .

Although the presence of hydrodynamic flow induces nonequilibrium conditions, the fluid has typically enough time to reach conditions of local equilibrium, and therefore a discussion based on free energies is admissible. We find that the total free energy can be expressed as

where the cavitation volume , α, β and γ are coefficients with dimensions of energy (see Methods section for details), τ the surface tension of the liquid–vapour interface, S N is the nematic order parameter in the absence of cavitation and ΔP is the absolute value of the pressure drop. Cavitation will occur when the free energy in equation (1) develops a minimum v*>0. This critical point is found from the two conditions ∂F tot /∂v=0 and ∂2F tot /∂v2=0 that determine the critical cavity volume and critical pressure

where is a reference pressure, , . For λ>0 and μ<0 a finite minimum develops.

The influence of the order parameter on cavitation can be found from studying how v* depends on S N ,

where . Thus, as S N grows, v* decreases and cavitation grows easier. This fact is confirmed by our computer simulations.