## Abstract

Pure liquids in thermodynamic equilibrium are structurally homogeneous. In liquid crystals, flow and light pulses are used to create reconfigurable domains with polar order. Moreover, through careful engineering of concerted microfluidic flows and localized optothermal fields, it is possible to achieve complete control over the nucleation, growth, and shape of such domains. Experiments, theory, and simulations indicate that the resulting structures can be stabilized indefinitely, provided the liquids are maintained in a controlled nonequilibrium state. The resulting sculpted liquids could find applications in microfluidic devices for selective encapsulation of solutes and particles into optically active compartments that interact with external stimuli.

## INTRODUCTION

Solid materials can simultaneously exhibit distinct structural phases, which can be manipulated to engineer functionality (*1*–*3*). Such structural phases, and the corresponding grain boundaries and defects, do not arise in pure liquids at equilibrium (*4*). On the one hand, liquids exhibit a number of attractive features, including short relaxation times, high diffusion coefficients, absolute compliance, and, of course, the ability to wet surfaces (*5*). On the other hand, however, imbuing pure liquids with additional functionality is challenging because of their inherent homogeneity. Complex behavior is generally encountered in multicomponent mixtures, be they synthetic or biological. In the particular case of biological systems, examples of self-organized, transient, and reconfigurable assemblies include raft domains, droplets, and other membraneless compartments (*6*). Such structures, however, are difficult to manipulate because they occur in out-of-equilibrium situations and generally involve multiple components that exhibit sharp miscibility gradients (leading to hydrophilic and hydrophobic domains). Recent efforts to “print” hydrophobic and hydrophilic domains into liquid mixtures by relying on the use of surfactant nanoparticles provide a notable demonstration of the possibilities afforded by gaining control over structure in nonequilibrium systems (*7*). Similarly, active matter with intrinsic topological defects is of particular relevance in this regard, where emergent structures can be created in the form of living colonies, tissues, and their biologically inspired synthetic counterparts (*8*–*12*). In such examples, intrinsic activity leads to motion and transitions between different rheological regimes.

Recent advances in the ability to control the competing effects of confinement and external fields by purposely designed micromanipulation tools have enabled seminal studies of nucleation, stability, and motion of topologically protected configurations in complex fluids (*13*–*17*), which can replace the need for multicomponent mixtures by creating distinct structural domains within a pure liquid. Liquid crystals (LCs) represent ideal systems for the study of spontaneous symmetry breaking, topological defects, orientational ordering, and phase transitions induced by applied stimuli (*18*). Even the simplest nematic phase, where the average orientation of rod-like molecules is characterized by a nematic director **n**, exhibits a wide range of switching mechanisms between uniformly aligned states (*19*). Nematic LCs (NLCs) can nucleate point and line defects (*20*) by rapid pressure or temperature quenches (*21*) and in the presence of colloidal inclusions (*22*–*25*). Recent work has shown that line defects can serve as microreactors in which to conduct polymerization reactions (*26*), offering intriguing prospects for future applications. Thermodynamic and anchoring transitions, textures, hydrodynamics, and flow instabilities in nematic mesophases have all been studied over the past decades. However, little is known about the coexistence and stability of driven orientational phases and the corresponding defects under geometric confinement. Recent results on nematic flows in microfluidic environments (*27*–*32*) have raised the possibility of tuning multistable defect patterns, and the transitions that arise between flow regimes, by controlling the shape of the channels, the anchoring conditions on channel walls, the temperature-dependent material properties, and, most importantly, the shear rate induced by the pressure differential that drives the flow.

## RESULTS

We demonstrate the creation and dynamic manipulation of defects and reconfigurable states in a pure NLC by simultaneous application of multiple external fields. Oriented polar-phase domains are generated and controlled through combinations of confinement, flow, and laser pulses. In contrast to previous reports, where only bulk states of both homeotropic and flow-aligned phases were considered, we are able to observe the critical behavior of the phase interface for the first time. Precision control over the production, growth rate, and shape of the domains is achieved through a combination of simulations and experiments that involve flow control and optothermal manipulation driven by laser tweezers. More specifically, a quantitative analytical model of the phase interface is proposed to predict the behavior of both phases from the geometric and material parameters of the experimental setup. The model is supported by detailed numerical simulations that include all relevant structural and hydrodynamic details. The insights provided by our theoretical and computational work are used to design a responsive system in which a metastable flow-aligned phase is repeatedly reconfigured by switching the flow direction, thereby permitting detailed study of the nucleation of solitons, domain walls, and point defects, followed by their relaxation dynamics as they seek to return to equilibrium.

Distinct domains of a flow-aligned phase of pentylcyanobiphenyl (5CB) NLC in a linear microchannel are nucleated through a temperature quench with a laser beam. As shown in Fig. 1A, the channel has a rectangular cross section, and its walls confer perpendicular surface alignment of the NLC (“homeotropic anchoring”) (see Materials and Methods). The director profile of the initial stationary state corresponds to a uniformly aligned homeotropic configuration along the *z* axis that appears black when observed between crossed polarizers. When the flow is turned on, the director remains predominantly aligned perpendicular to the substrate but is slightly deflected in the flow direction due to flow alignment, changing the birefringent appearance from black to bright colors that depend on the flow velocity (Fig. 1B). We name this flow regime the “bowser state,” after the bowed shape of the director profile in contrast with the flow-aligned state, which we will refer to as the “dowser state.” A black isotropic island is created where the NLC is heated into the isotropic phase by the laser tweezers; after the light is switched off and the NLC is quenched into the nematic phase, the initial tangle of defects relaxes into a flow-aligned state, bounded by a disclination loop (Fig. 1B). The flow-aligned domains then evolve with flow (Fig. 1C) and can either grow or annihilate depending on the flow velocity (movies S1 and S2). The flow-oriented phase, whose director field in the midplane of the channel is aligned horizontally, changes direction by a half-turn from the bottom to the top of the channel. This system represents an ideal experimental model of a quasi–two-dimensional (2D) orientational phase, described by a unit vector field of the midplane director alignment, as recently reported in the strong flow regime (*v* > 100 μm/s) in nematodynamics (*30*), and it is analogous to the so-called dowser field in nematostatics, after which we adopted the name (*33*).

The dynamics of the director field, described by the in-plane angle ϕ relative to the *x* axis, is driven by the elasticity and the flow alignment. Assuming uniform flow with midplane velocity *v* along the *x* axis and neglecting splay-bend anisotropy, we obtain a linear coupling of director and flow velocity(1)

This coupling is reminiscent of that of the dowser field to the thickness gradient (*33*–*35*). The equation can be recognized as a damped sine-Gordon equation with characteristic length , where *K* is the single elastic constant, γ_{1} and γ_{2} are viscosity parameters, and *h* is the height of the channel. Details on the general derivation, including the velocity gradient term, which is essential for describing flow in channel junctions (*32*), are provided in Materials and Methods.

The flow-aligned dowser state is stable under strong flows but unstable in weak flows. The dowser domains are seen to grow and shrink, depending on the flow velocity, as shown both in experiments (Fig. 2, A, C and F, and movie S2) and in numerical simulations (Fig. 2, B, D, and E, and movie S3). We have constructed a quantitative phenomenological model for the dynamics of the domain size (see Materials and Methods). In such a model, the growth rate of a circular domain with radius *r* is given by(2)where *T* is the disclination line tension and γ_{r} is the viscosity parameter associated with the drag force on moving disclination lines, which are both assumed to be constant. The critical velocity is . The stability of a flow-aligned (ϕ = 0) dowser domain is conditioned by the midplane velocity in the channel, as summarized in the phase diagram of Fig. 2H. Below *v*_{c}, all dowser domains are unstable and gradually shrink in size until they disappear. Above *v*_{c}, the behavior of dowser domains depends on their initial size. If the domain size is smaller than the critical radius , the domain shrinks over time and annihilates. If *r* > *r*_{c}, the domain grows until it reaches the side walls of the channels and expands along the channel from then on (Fig. 2, A and B).

Equation 2 can be integrated analytically to yield the time dependence of the flow-aligned dowser domain size(3)where and *r*_{0} is the initial radius of the loop. We have fitted Eq. 3 to the experimental data in Fig. 2F and fig. S1 through parameters *a*, *r*_{0}, and *r*_{c} and obtained good agreement with the model. From the fitting parameters, a *v*_{c} of (56.4 ± 1.4) μm/s was determined as the point where the inverse of the *r*_{c} reaches zero (Fig. 2G). Fitting over parameter *a* yields a similar value (56.8 ± 1.2) μm/s (fig. S2). The critical velocity, calculated directly from the dimensions of a channel and the viscoelastic properties of 5CB, is ≈ 42 μm/s. The agreement with both values obtained from the fit is reasonable, particularly considering the simplifying assumptions involved in the theoretical model. A separate analysis, performed on channels of different sizes, reveals a comparable behavior and fits the model reasonably well (fig. S3).

The dowser domains can be manipulated in various ways through careful application of the laser tweezers. A preexisting bulk dowser state may be created upstream by prior quenching or simply by the initial conditions of the nematic at the influx. One can produce a steady stream of domains by dissecting the original bulk dowser with a moving laser spot (Fig. 3A and movie S4) whose role is to constantly melt the sides of a phase boundary. A growing domain at higher flow velocity can be longitudinally split in half by a static laser beam at lower light intensities (Fig. 3B and movie S5). One can observe changing birefringent colors as the domain traverses a light-generated obstacle. The laser tweezers therefore enable dynamic control over the size, number, and lifetime of generated dowser domains, which can be further manipulated by periodic flow velocity modulations. In a uniform flow, the dowser field aligns uniformly along the flow direction and either grows or shrinks, depending on the velocity regime. By careful tuning and active control of the flow, a constant-size domain can be maintained over a period of tens of seconds (Fig. 3C). Here, the flowing dowser domain is stabilized in an oscillatory flow, maintained by a sinusoidal modulation of the driving pressure around the target flow rate (70 μm/s < *v* < 95 μm/s, with a period of 3.5 s).

Our model, as expressed in Eq. 1, predicts that, when the flow direction is reversed, the previous equilibrium state becomes the least favorable, leading to a rapid reversal of orientation. Figure 3D and movie S6 show a dowser domain in an alternating flow where flow reversal can be observed twice per period (−90 μm/s < *v* < +90 μm/s with a period of 2.0 s). When the velocity drops below the critical value, the domain begins to shrink. As the flow reverses, the reorientation begins at the leading and trailing points of the domain, creating reorientation fronts that propagate inward (Fig. 3D). The part with the old orientation is energetically unfavorable, causing it to shrink to a narrow 2π soliton. When the velocity passes *v*_{c}, the domain starts to grow again, leading to disappearance of solitons through discontinuous director rearrangement (fig. S4). The soliton size is set by the characteristic length ξ from the sine-Gordon equation, which becomes very small under strong flows.

The line tension equation Eq. 2 holds locally and describes the curvature flow of the dowser domain boundary. When the dowser orientation ϕ is not aligned with the flow, the domain growth rate varies along the boundary and may even change sign. Figure 3D shows invaginations that appear in the dowser domain, where the soliton extends to the boundary. The boundary of the domain is thus controlled by the internal orientation, leading to noncircular shapes that are reminiscent of living cells during division. In Fig. 3E, one can appreciate that bulbous parts appear on the shrinking domains; these are due to small residual flows that cause a difference in line tension between regions with different director orientations.

The dowser state, being a polar unit vector field, supports the existence of point defects with an integer winding number, which appear in pairs, connected by a soliton. As shown in Fig. 3E and movie S7, turning off the flow increases the characteristic length with decreasing velocity, allowing us to observe the detailed structure of the soliton. In this simple case, the sine-Gordon equation can be solved analytically, giving the transverse profile of the soliton as ϕ(*y*) = 4 arctan *e*^{±y/ξ} (inset in Fig. 3E), analogous to the soliton profile seen in the dowser field under the cuneitropism effect (*33*). In addition to defect pairs, closed circular solitons may appear during oscillatory flow (fig. S4).

## DISCUSSION

In close analogy to other heterogeneous mixtures with phase boundary interfaces, such as fogs and aerosols, microfluidic droplets, fluids at the triple point, and the nematic-isotropic phase transition itself, surface tension plays a crucial role in controlling the “evaporation” and growth of the low-energy phase, depending on the curvature of the interface. In the system presented here, the phases have different symmetries, are topologically incompatible, and exhibit what is essentially a first-order phase transition. The dowser state is orientationally anisotropic, with its own elastic behavior, topological defects, and solitons, whereas the bowser state is effectively isotropic and simple in the simplified 2D view. The dowser-phase orientation couples with the velocity field and can have either lower or higher energy than the bowser, depending on its orientation. In this way, the shape, splitting, and coalescence of domains can be controlled.

As the dowser field also couples to external magnetic and electric fields, as well as the gradients of the channel thickness, we suggest that combinations of these external cues could offer exceptional possibilities for shape control, flow steering, and optical tuning. As the response to the external stimuli is directly and clearly observable through birefringence, it could also be used to measure the viscoelastic and rheological properties of the material itself. The dowser domains, readily produced in large quantities and sizes with laser tweezers, are analogous to droplets and shells and can serve similar purposes. It should be possible, for example, to create nested domains for the production of 2D shells with an enclosed volume of the dowser phase. It should also be possible to conduct chemical reactions in such domains or sequester different components according to orientational or dielectric affinity (as opposed to hydrophobic/hydrophilic contrast), as recently shown for defects around microparticles (*26*). One could envision a 3D printing system for liquids in which complex out-of-equilibrium structures are created and stabilized by relying on the principles outlined here. Last, the models obtained here from experiments with standard thermotropic LCs could be readily applied to active and biological materials with nematic behavior. Controlling alignment and separation into isolated orientational domains may be seen as a technical tool with potential for applications in biophysics, chemistry, and chemical engineering.

## MATERIALS AND METHODS

### Materials and experimental procedures

We used a single-component nematic material 5CB (SYNTHON Chemicals), which has the nematic phase in temperature range 18°C < *T* < 35°C, in all of our experiments. The microfluidic channels had a rectangular cross section, with height *h* ≈ 12 μm, width *w* = 100 μm, and length *L* = 20 mm. The channels were fabricated out of polydimethylsiloxane (PDMS; SYLGARD 184, Dow Corning) reliefs and indium tin oxide (ITO)–coated glass substrates (Xinyan Technology) by following standard soft lithography procedures (*30*). The ITO coating was used as an absorber of the infrared (IR) laser light and provided very good control for the local heating of 5CB. The channel walls were chemically treated with 0.2 weight % aqueous solution of silane dimethyl-octadecyl-3-aminopropyl-trimethoxysilyl chloride (abcr GmbH) to induce strong homeotropic surface anchoring for 5CB molecules (*30*). The microfluidic channels were filled with 5CB in its isotropic phase and allowed to cool down to nematic phase at room temperature before starting the flow experiments. We drove and precisely controlled the fluid flow by using a pressure-driven microfluidic flow control system (OB1, Elveflow). We applied and varied flow rates in the range 0.05 to 1.50 μl/hour corresponding to a flow velocity *v*, ranging from 5 to 250 μm/s. In some experiments, the nematodynamics was controlled by adjusting time-dependent flow driving using built-in flow profile routines of the microfluidic controller. The characteristic Reynolds number *Re* = ρ*vl*/η for 5CB having effective dynamic viscosity (*32*) η ≈ 50 mPa⋅s ranged between 10^{−6} and 10^{−4}. Here, ρ ≈ 1.024 kg/m^{3} is the material’s density, and *l* = 4*wd*/2(*w* + *h*) ≈ 21 μm is the hydraulic diameter of the channels. The corresponding Ericksen number *Er* = η*vl*/*K*, with *K* = 5.5 pN being the 5CB single elastic constant approximation, varied between 0.8 and 40. All the experiments were conducted at room temperature.

### Polarized light microscopy, laser tweezers, and image acquisition

The flow regimes, reorientation dynamics, and flow-driven deformations of 5CB in microchannels with homeotropic surface anchoring were studied by polarized light microscopy (Nikon, Eclipse Ti-U, equipped with CFI Plan 2× and 10× objectives). The samples were observed between crossed polarizers in transmission mode. In addition, we used a laser tweezers setup build around the inverted optical microscope with an IR fiber laser operating at 1064 nm as a light source and a pair of acousto-optic deflectors driven by a computerized system (Aresis, Tweez 200si) for precise laser beam manipulation. The laser power was varied between 20 and 200 mW in the sample plane, and the Gaussian beam profile was primarily used for local heating of the NLC above its clearing temperature. Full high-definition color videos were recorded at a frame rate of 30 frames per second using a digital complementary metal-oxide semiconductor camera (Canon, EOS 750D), attached by a C-mount compatible adapter (LMScope) to the microscope. The image analysis was performed using the software ImageJ.

### Numerical simulation details

The bulk free energy of the NLC, *F*, is defined as(4)where *f*_{LdG} is the short-range free energy, *f*_{el} is the long-range elastic energy, *f*_{D} is the laser-induced dielectric interaction energy, and *f*_{surf} is the surface free energy due to anchoring. *f*_{LdG} is given by a Landau–de Gennes expression of the form (*18*)(5)

Parameter *w* controls the magnitude of *q*_{0}, namely, the equilibrium scalar order parameter via . The elastic energy *f*_{el} is written as (using Einstein summation rule)(6)where ∂_{i} indicates a spatial derivative over the *i*th coordinate. If the system is uniaxial, the *L*’s in Eq. 6 can be mapped to the Frank elastic constants *K*’s via(7)

By assuming a single elastic constant *K*_{11} = *K*_{22} = *K*_{33} = *K*_{24} ≡ *K*, one has and *L*_{2} = *L*_{3} = *L*_{4} = 0. Pointwise, **n** is the eigenvector associated with the greatest eigenvalue of the **Q**-tensor at each lattice point. The free energy associated with anisotropic dielectric constants readswhere ε_{0} is the vacuum permittivity constant and ε_{ij} is the dielectric permittivity tensor related to the **Q**-tensor asin which ε_{||} and ε_{⊥} are the permittivities parallel and perpendicular to the nematic director, respectively.

To simulate NLC’s flowing dynamics, a hybrid lattice Boltzmann method is used to simultaneously solve a Beris-Edwards equation and a momentum equation that accounts for the backflow effects. By introducing a velocity gradient *W*_{ij} = ∂_{i}*v*_{j}, a strain rate , a vorticity tensor , and a generalized advection term(8)one can write the Beris-Edwards equation (*36*) according to(9)where ∂_{t} is a partial derivative over time. The constant ξ is related to the material’s aspect ratio and relates to the alignment parameter in the Ericksen-Leslie-Parodi theory . Γ is related to the rotational viscosity γ_{1} of the system by (*37*). The molecular field **H**, which drives the system toward thermodynamic equilibrium, is given by(10)where […]^{st} is a symmetric and traceless operator. When velocity is absent, i.e., **v**(**r**)≡0, the Beris-Edwards equation (Eq. 9) reduces to the Ginzburg-Landau equation

To calculate the static structures of ±1/2 defects, we adopted the above equation to solve for the **Q**-tensor at equilibrium.

Homeotropic anchoring was implemented through a Rapini-Papoular expression (*18*) that penalizes any deviation of the **Q**-tensor from **Q**_{0}, namely, a surface-preferred **Q**-tensor. The associated free-energy expression is given by(11)

The evolution of the surface **Q**-field is governed by (*31*)(12)where Γ_{s} = Γ/ξ_{N} with , namely, nematic coherence length. The above equation is equivalent to the mixed boundary condition given in (*38*) for steady flows.

The momentum equation for the nematics is written as (*37*)(13)

The stress Π is defined as(14)where η is the isotropic viscosity, and the hydrostatic pressure *P*_{0} is given by(15)

The temperature *T* is related to the speed of sound *c*_{s} by . We solve the evolution equation (Eq. 9) using a finite difference method. The momentum equation (Eq. 13) is solved simultaneously via a lattice Boltzmann method over a D3Q15 grid (*39*). The implementation of stress follows the approach proposed by Guo *et al*. (*40*). Our model and implementation were validated by comparing our simulation results to predictions using the Ericksen-Leslie-Parodi theory (*41*). The units are chosen as follows: The unit length *a* is chosen such that unit length *a* = ξ_{N} = 7 nm, viscosity γ_{1} = 0.07 Pa⋅s, and elastic constants splay *K*_{11} = 6 pN, twist *K*_{22} = 3.9 pN, bend *K*_{33} = 8.2 pN, and saddle-splay *K*_{24} = 7 pN, mimicking the material properties of 5CB. We refer the reader to (*31*) for additional details on the numerical methods used here.

### Dowser field orientation

To derive an effective 2D theory of a dowser state, we took the ansatz for the dowser profile of the director field(16)where the unit vector **d** is a dowser orientational field in the *xy* plane, *z* is the vertical position in the channel of height *h*, and **e**_{z} is a unit vector along the *z* direction. A Poiseuille flow profile is assumed(17)where **v** is a vector field in the *xy* plane of the channel. Given that we are interested only in the behavior of the director field, we write the dissipation function, omitting the terms that do not include time derivatives of **n** (*42*)(18)where and are viscosity coefficients in the Ericksen-Leslie-Parodi formulation of nematodynamics and a dot indicates the time derivative, assuming that entire vertical profile advects with the midplane velocity **v**. Total dissipation in the system is given as a volume integral over the dissipation function *D*. Since the director and the velocity profile in the *z* direction are fixed by Eqs. 16 and 17, respectively, we are free to perform the integration over the *z* axis, obtaining effectively a dissipation function *D*_{2D} = ∫*D*d*z* of 2D processes(19)

Next, we write the elastic free-energy density of a dowser structure in a single elastic constant (*K*) approximation , which can again be integrated over the *z* axis *f*_{2D} = ∫*f*_{el}d*z*, obtaining(20)

We write dowser field as **d** = (cos ϕ, sin ϕ) and follow the Lagrange formalism for ϕ angle , obtaining the master equation for the dowser orientation(21)

### Effective free energy of dowser domains

The time derivative of the dowser orientation can also be viewed as a relaxation under the effective potential *U* that depends on the velocity field. The equation of motion in that case is(22)

The master equation for is recovered for(23)

Integrating the above equation over ϕ leads to(24)which can be written in a covariant form(25)where a constant *C* is not dependent on *ϕ*. *C* has to vanish at **v** = 0 and must be at least quadratic in *v* to preserve the invariant form. Since we are interested only in the linear response to the velocity field, we can set *C* = 0.

Using the effective potential *U*, we can phenomenologically construct an effective free energy *F* of a dowser state in microfluidic confinement in contact with a homeotropic nematic state (with ansatz **n** = **e**_{z})(26)where *T* is the line tension of a nematic disclination (*18*). Specifically, we are interested in the free energy of circular dowser domains with homogeneous alignment of the dowser vector **d** in a flow field that is homogeneous in the *xy* plane. This substantially simplifies the expression for the free energy(27)where . The dynamics of the loop growth or annihilation is given by(28)where the viscosity parameter γ_{r} is due to a drag force on a moving disclination line. The line tension contribution becomes dominant at small radii, leading to universal annihilation behavior of shrinking loops (fig. S1).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/2/eaav4283/DC1

Fig. S1. Universality of the dowser domain annihilation dynamics at small radii.

Fig. S2. Velocity dependence of parameter *a* from Eq. 3 for the dowser domain size, as extracted from Fig. 2F.

Fig. S3. An independent study, performed in 400-μm-wide and 15-μm-deep channels, shows equivalent behavior of dowser domain dynamics as is discussed in the main text.

Fig. S4. Dowser field relaxation dynamics.

Movie S1. Laser tweezers-induced nucleation of the dowser domains in a pressure-driven nematic microflow.

Movie S2. Expansion and contraction of laser-nucleated dowser domains in a moderate nematic microflow.

Movie S3. Growing and shrinking dowser domains in numerically simulated nematic microflows.

Movie S4. A steady stream of dowser domains is produced by chopping the bulk dowser state with a moving laser spot.

Movie S5. A growing dowser domain is longitudinally split into two by a static laser spot.

Movie S6. Dowser domain reconfiguration under an oscillatory flow.

Movie S7. Relaxation dynamics of dowser domains after shutting off the flow.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank J. A. Martinez-Gonzalez and M. Ravnik for stimulating discussions on topics related to this work.

**Funding:**This work was conducted with support from the Slovenian Research Agency (ARRS) under contracts P1-0055 (to T.E. and U.T.), P1-0099 (to Ž.K. and S.Č.), P1-0192 (to N.O.), J1-6723 (to U.T.), L1-8135 (to Ž.K., S.Č., and U.T.), and J1-9149 (to S.Č.), and with support from the National Science Foundation of the United States of America under grant DMR 1710318 to J.J.d.P. U.T. would like to thank COST Action MP1205 “Advances in Optofluidics: Integration of Optical Control and Photonics with Microfluidics” and COST Action MP1305 “Flowing Matter” for supporting his research activities. R.Z. is grateful for the support of the University of Chicago Research Computing Center for assistance with the simulations carried out in this work.

**Author contributions:**U.T. and J.J.d.P. designed the research. T.E. and U.T. conducted the experiments. R.Z. and J.J.d.P. performed the numerical simulations and analyzed the data. Ž.K. and S.Č. developed the theoretical model and analyzed the results. N.O. contributed measurements for the Supplementary Materials. U.T., S.Č., Ž.K., R.Z., and J.J.d.P. wrote the paper. U.T. and J.J.d.P. supervised the research. All authors discussed the progress of the research and contributed to the final version of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).