Bifurcation dynamics of a particle-encapsulating droplet in shear flow

TL;DR

This paper investigates the complex bifurcation behavior of a particle encapsulated within a deforming droplet under shear flow, revealing new equilibrium states and the underlying hydrodynamic interactions through a dynamical systems approach.

Contribution

It introduces a detailed numerical analysis of bifurcations in particle-encapsulating droplets, highlighting new equilibrium solutions and the role of hydrodynamic forces in their emergence.

Findings

Identification of multiple equilibrium solutions including eccentric and time-periodic states.

Discovery of bifurcation points such as supercritical pitchfork and Hopf bifurcations.

Elucidation of hydrodynamic interactions driving the bifurcations.

Abstract

To understand the behavior of composite fluid particles such as nucleated cells and double-emulsions in flow, we study a finite-size particle encapsulated in a deforming droplet under shear flow as a model system. In addition to its concentric particle-droplet configuration, we numerically explore other eccentric and time-periodic equilibrium solutions, which emerge spontaneously via supercritical pitchfork and Hopf bifurcations. We present the loci of these solutions around the codimenstion-two point. We adopt a dynamical system approach to model and characterize the coupled behavior of the two bifurcations. By exploring the flow fields and hydrodynamic forces in detail, we identify the role of hydrodynamic particle-droplet interaction which gives rise to these bifurcations.

Equations

$${\dot{\rho }}_1$$

$${\dot{\rho }}_2$$

$${\dot{\xi }}_1$$

$${\dot{\xi }}_2$$

$${\mu }_1$$

$${\mu }_2$$

Full text

††preprint:

Droplets, capsules and vesicles in flow often exhibit interestingly rich dynamics even in the linear shear flow Stone (1994); Smith et al. (2004); Sibillo et al. (2006); Skotheim and Secomb (2007); Omori et al. (2012); Kraus et al. (1996); Misbah (2006); Noguchi and Gompper (2007); Deschamps et al. (2009). Despite the substantial work on the dynamics of these soft systems enclosing homogeneous fluids, limited effort has been directed to studying their behavior when they include an internal structure. However, such a configuration is common in nature and engineering applications: cells like leukocytes, and megakaryocytes contain nucleus up to $50\sim 80\%$ of themselves in volume Turgeon (2005); double-emulsions playing an important role in chemical and pharmaceutical engineering are featured with a core-shell geometry Stone and Leal (1990); Utada et al. (2005); Guzowski et al. (2013) ; droplet-based encapsulation for high-throughput biological assays utilizes droplets as micro-chambers to compartment cells for analysis at the single-cell level, where the cell size can be comparable to the droplet size in certain applications  He et al. (2005); Chabert and Viovy (2008); Mazutis et al. (2013).

These fluid particles are characterized by complex hydrodynamic interactions between the internal structures and the external interface. Few works conducted for nucleated model cells in shear  Veerapaneni et al. (2011); Kaoui et al. (2013); Levant and Steinberg (2014); Luo et al. (2015) all assumed their compound structures to be concentric, preserving the rotational symmetry of order $2$ (C$2$) about the $y$ axis and reflection symmetry about the $y=0$ shear plane (see Fig. 1a). The symmetries do hold for a single shear-driven deformable particle which attains a steady ellipsoidal shape undergoing tank-treading motion Rallison (1984); Seifert (1997); Barthès-Biesel (2016). Yet, they are not guaranteed in the presence of an internal structure.

In this Letter we focus on the stability of the concentricity of composite fluid particles. By considering a droplet encaging a spherical particle as a model system, we formulate the following questions: will the composite structures remain concentric? How does the dynamics depend on interfacial tension and particle size? What is the role of the hydrodynamic interaction?

We begin our discussion by presenting $3$D hydrodynamic simulations of a compound particle-droplet subjected to unbounded shear $\mathbf{U}^{\infty}=\mathbf{G}\cdot\mathbf{x}$, in the creeping flow regime, where the only non-zero component $G_{xz}=\dot{\gamma}$ represents the shear rate (Fig. 1a). The incompressible Stokes equations are solved by a boundary integral method (see Supplemental material and  Reigh et al. (2017) for details). The immiscible Newtonian fluids inside and outside the droplet have the same viscosity $\eta$; its surfactant-free interface has a uniform surface tension $\sigma$. The particle has a no-slip surface, freely translating and rotating subject to zero hydrodynamic force and torque. The droplet interface satisfies the standard stress balance condition bc_ ; Leal (2007). The radii of the particle and the undeformed droplet are $a$ and $R$ respectively; the size ratio is denoted by $\alpha=a/R$ with $\alpha\in\left(0,1\right)$. The capillary number $\mathrm{Ca}=\eta\dot{\gamma}R/\sigma$ indicates the ratio between viscous forces and capillary forces, limited to the regime without droplet break-up. All length scales are scaled by $R$.

We initially displace the particle away from the droplet center by a perturbative offset $\mathbf{d}=\left(d_{x},d_{y},d_{z}\right)$, then focusing on the time evolution of $d_{y}$ and $d_{xz}=\sqrt{d_{x}^{2}+d_{z}^{2}}$ representing the spanwise and in-plane displacements respectively. $\Delta_{y}$ and $\Delta_{xz}$ denote their equilibrium values when the system reaches a steady or time-periodic state.

We show the evolution of displacements in Fig. 2, presenting three typical $\mathrm{Ca}$-dependent scenarios for a particle of size ratio $\alpha=0.4$ (see Supplemental videos). When $\mathrm{Ca}=0.2$, the in-plane displacement $d_{xz}$ decays asymptotically to zero after a transient growth while the spanwise offset $d_{y}$ increases to a saturated value $\Delta_{y}\approx 0.48$ indicating the broken reflection symmetry. The particle rotates steadily near the lateral edge of the droplet interface (Fig. 2b). Increasing $\mathrm{Ca}$ to $0.275$, $\Delta_{y}$ decreases to $0.21$ approximately, while $d_{xz}$ reaches a time-periodic equilibrium cycle with a maximum of $\Delta^{\mathrm{max}}_{xz}\approx 0.11$ and a minimum of $\Delta^{\mathrm{min}}_{xz}\approx 0.04$. The particle follows an orbital trajectory on the $y=\Delta_{y}$ plane as it reaches a limit cycle solution in the $\left(x,z\right)$ space (Fig. 2d, f), implying that the C$2$ symmetry and time invariance are also broken. At $\mathrm{Ca}=0.29$, the system recovers steadiness and concentricity, $\Delta_{y}=\Delta_{xz}=0$. These scenarios suggest the appearance of bifurcating solutions by reducing $\mathrm{Ca}$: above a critical value $\mathrm{Ca}_{\mathrm{c}}\left(\alpha\right)$, the composite system stays concentric, corresponding to a stable fixed point solution; it bifurcates across $\mathrm{Ca}_{\mathrm{c}}\left(\alpha\right)$ towards a steady spanwise migration (SM) and/or in-plane orbiting (IPO) motion.

An investigation spanning the $\left(\mathrm{Ca},\alpha\right)$ space further reveals that these two modes, i.e. SM and IPO, appear spontaneously through supercritical pitchfork and Hopf bifurcations respectively. To study the evolution of the two modes individually, we perform decoupled simulations with kinematic constraints of either $d_{xz}=0$ (pure SM) or $d_{y}=0$ (pure IPO). Their corresponding equilibrium displacements $\tilde{\Delta}_{y}$ and $\tilde{\Delta}^{\mathrm{min}/\mathrm{max}}_{xz}$ are shown in Fig. 3. For all $\alpha$, we observe $\tilde{\Delta}_{y}$ decreases with $\mathrm{Ca}$, becoming zero when $\mathrm{Ca}$ exceeds a critical value $\mathrm{Ca}_{\mathrm{c}}\left(\alpha\right)$, so does $\tilde{\Delta}^{\mathrm{min}/\mathrm{max}}_{xz}$. They both vary quadratically in the vicinity of their corresponding $\mathrm{Ca}_{\mathrm{c}}\left(\alpha\right)$. This is confirmed by the linear fitting of $\tilde{\Delta}_{y}$ and, $\tilde{\mathcal{A}}_{xz}=\tilde{\Delta}^{\mathrm{max}}_{xz}-\tilde{\Delta}^{\mathrm{min}}_{xz}$ indicating the oscillating amplitude, versus $\left[\mathrm{Ca}_{\mathrm{c}}\left(\alpha\right)-\mathrm{Ca}\right]^{1/2}$ (Fig. 3d), where $\mathrm{Ca}_{\mathrm{c}}\left(\alpha\right)$ is obtained simultaneously. The successful fitting passing through the origin verifies the emergence of the two bifurcations. It is worth-pointing that broken reflection symmetry by SM is indeed the signature of pitchfork bifurcation, so as broken time invariance by IPO of Hopf bifurcation.

Let us return to the constraint-free cases, where the nonlinear interaction of the two modes results in a more complex dependence of the equilibrium solutions on $\mathrm{Ca}$ and $\alpha$ (Fig. 3). For $\alpha=0.3$, the in-plane amplitudes (cyan) reach their maxima around $\mathrm{Ca}=0.3$, from where they decrease almost linearly/quadratically with decreasing/increasing $\mathrm{Ca}$. In their quadratic parts, the coupled (cyan) and decoupled amplitudes (green) overlap in the vicinity of their common critical point. In contrast, the spanwise offsets (purple) are larger than those of the decoupled cases (red). For $\alpha=0.4$, the spanwise offsets coincide precisely with the decoupled counterparts for all $\mathrm{Ca}$, while the in-plane amplitudes exhibit non-monotonic $\mathrm{Ca}$-dependence as for $\alpha=0.3$ and they are below the decoupled values. For $\alpha=0.5$, perfect coincidence between the spanwise offsets also holds as in the $\alpha=0.4$ case, while $\Delta^{\mathrm{min}/\mathrm{max}}_{xz}\equiv 0$ for all $\mathrm{Ca}$, i.e., the Hopf bifurcation is inhibited.

The complexity is better unraveled by the parametric portrait quartering the $\left(\mathrm{Ca},\alpha\right)$ parameter space into the following solution types (Fig. 4): ’concentric’ implying the absence of both modes, ’pure IPO’, ’pure SM’, and ’mixed’ indicating the coexistence of both modes. A codimension-two point $\left(\mathrm{Ca}_{\ast},\alpha_{\ast}\right)\approx\left(0.286,0.4\right)$ is pinpointed at the intersection of the two marginal curves $H_{1}^{\pm}$ (circle) and $P_{2}^{\pm}$ (triangle) which correspond to the Hopf and pitchfork bifurcations, respectively. Other branches bifurcating from this point are $T_{1}$ separating ’pure SM’ and ’mixed’, $T_{2}$ separating ’pure IPO’ and ’mixed’. Note + and - denote the upper and lower branches of the marginal curves.

We now interpret the bifurcation in the neighborhood of the codimension-two point $\left(\mathrm{Ca}_{\ast},\alpha_{\ast}\right)$ based on a normal-form analysis. By coupling the amplitude equations of the Hopf and pitchfork bifurcations, we obtain a normal form similar to that of the Hopf-Hopf bifurcation in Ref. Kuznetsov (2013) where the amplitudes are independent of phase evolution. Denoting the square of the in-plane and spanwise amplitudes by $\rho_{1}=\left(\Delta^{\mathrm{max}}_{xz}-\Delta^{\mathrm{min}}_{xz}\right)^{2}$ and $\rho_{2}=\Delta_{y}^{2}$, the truncated amplitude system is expressed as

[TABLE]

where $\mu_{i}$ represent the linear growth-rates of the individual modes and $p_{ij}$ the nonlinear coupling coefficients. Because of the supercritical nature of the two bifurcations, $p_{11}<0$ and $p_{22}<0$. Physically, the amplitudes tend to asymptotic values with decreasing $\mathrm{Ca}$ (see Fig. 3) owing to the confinement of droplet. Because $p_{11}p_{22}>0$, our problem is in the category of the so-called simple cases, for which we have neglected $\rho_{1}\rho^{2}_{2}$ and $\rho_{2}\rho^{2}_{1}$ without changing the bifurcation topology Kuznetsov (2013). By introducing new phase variables $\xi_{1}=-p_{11}\rho_{1}$ and $\xi_{2}=-p_{22}\rho_{2}$, we obtain

[TABLE]

where $\theta=p_{12}/p_{22}$ and $\delta=p_{21}/p_{11}$. Applying at leading order the affine transformation

[TABLE]

in the vicinity of $\left(\mathrm{Ca}_{\ast},\alpha_{\ast}\right)$, we map the parameter space from $\left(\mathrm{Ca},\alpha\right)$ to $\left(\mu_{1},\mu_{2}\right)$ (inset of Fig. 4b), where $k_{1}\approx-2.3$ and $k_{2}\approx-40$ denote the slope of $H_{1}$ and $P_{2}$ curves at $\left(\mathrm{Ca}_{\ast},\alpha_{\ast}\right)$, with $C\approx 0.043$ derived from the growth rate of $\rho_{i}$. The slopes of $T_{1}$ and $T_{2}$ further determine $\theta\approx 0.71$ and $\delta\approx-0.69$. The parametric portrait therefore corresponds to case III described in Ref. Kuznetsov (2013), characterized by six regions:
$\mathbf{1}$  corresponding to ’concentric’;
$\mathbf{2}$  to ’pure IPO’;
$\mathbf{3}$  and
$\mathbf{4}$  to ’mixed’ separated by $P_{2}^{+}$;
$\mathbf{5}$  and
$\mathbf{6}$  to ’pure SM’ separated by $H_{1}^{+}$ (see Supplemental material for their phase portraits).

The parametric portrait and normal form both reveal nonlinear mode interactions as a fingerprint of the present bifurcation. In the absence of SM, IPO appears stably in regions
$\mathbf{2}$  $\sim$
$\mathbf{5}$  but it is suppressed due to the nonlinear interaction with SM, as reflected by the phase portraits of
$\mathbf{3}$  $\sim$
$\mathbf{5}$  all including an unstable saddle-node equilibrium $\left(\mu_{1},0\right)$, as well as by the sign of $\theta$. Consequently, pure IPO only survives in
$\mathbf{2}$  . Besides, without IPO, pure SM is stable in regions
$\mathbf{4}$  $\sim$
$\mathbf{6}$  , while IPO promotes SM to expand its locus further to
$\mathbf{3}$  that indeed involves a stable equilibrium $\left(\mu_{1},\mu_{2}\right)$. This promotion results from the sign of $\delta$.

We next reveal the mechanisms underlying the bifurcations, firstly focusing on $H_{1}$ and $P_{2}$ separately. The shear flow can be decomposed into a rotational and extensional part, and we found that the former alone does not contribute to the particle’s cross-stream motions. The eccentricity is mostly driven by the extensional part $\mathbf{U}^{\infty}_{\mathrm{E}}\left(\mathbf{x}\right)=\mathbf{E}\cdot\mathbf{x}$ with $E_{xz}=E_{zx}=\dot{\gamma}/2$. For a system with imposed concentricity in $\mathbf{U}^{\infty}_{\mathrm{E}}\left(\mathbf{x}\right)$, Fig. 5 displays the droplet-induced disturbance flows $\hat{\mathbf{u}}_{\mathrm{d}}$ on $y=0$ (a) and $x=z$ (b) plane, which preserve reflection symmetries about each other. The major/minor axis of the ellipsoid-shaped droplet lies on $x=\pm z$ plane. The disturbance flow is induced to satisfy zero normal velocities on the interface. On the shear plane, it approaches/leaves the origin along the major/minor axis. It resembles the stagnation point flow, where the origin is kinematically unstable. This initiates the in-plane motion and indeed the particle moves along the minor axis and eventually touches the droplet for any $\mathrm{Ca}$ if we free its in-plane motion. This scenario is altered by the rotational flow, which relocates the particle between the two axes cyclically. Consequently, it is centralized/decentralized by the inward/outward flow after every relocation. The inward and outward flows roughly balance at $\mathrm{Ca}=0.1$; while the former dominates the latter at $\mathrm{Ca}=0.2$, hence overcoming the kinematic instability and leading to a concentric preference. This might explain the quenching of IPO when $\mathrm{Ca}$ increases across the marginal $H_{1}$ curve. On the $x=z$ plane (Fig. 5b), the flow resembles a parallel compressional flow which reaches the maximum strength at $y=0$ and weakens in $\pm y$ directions. When the particle undergoes a spanwise $d_{y}$ perturbation (say $d_{y}>0$), it experiences the strongest compression on its lower part and the $y$-gradient of that compressional flow will produces a viscous shear force in the spanwise direction (see discussion below) that further amplifies this perturbation, triggering the pitchfork bifurcation. Note that the flow on the $x=-z$ plane may conversely help centralize the particle, yet, it is weaker for any $\mathrm{Ca}>0$. Moreover, we conduct simulations fixing a certain spanwise offset $d_{y}^{\mathrm{fix}}>0$ with $d_{x}=d_{z}=0$ for the shear flow, recording the spanwise hydrodynamic forces $F_{Y}=F_{Y}^{p}+F_{Y}^{\eta}$ on the particle (Fig. 5c), where $F_{Y}^{p}$ (resp. $F_{Y}^{\eta}$) represents the pressure (resp. viscous) contribution which centralizes (resp. decentralizes) the particle. As shown, the viscous shear force $F_{Y}^{\eta}|_{s}$ accounts for the major contribution to $F_{Y}^{\eta}$, supporting the above arguments of compression-induced viscous destabilization. The pressure force becomes stronger with $\mathrm{Ca}$ and dominates the viscous part when $\mathrm{Ca}$ exceeds a critical value. In fact, the droplet with larger $\mathrm{Ca}$ displays a lateral protrusion accompanying a local curvature increase (Fig. 5e), generating a stronger pressure to center the particle. This clarifies why SM vanishes when $\mathrm{Ca}$ crosses $P_{2}$ curve. Upon having elucidated $H_{1}$ and $P_{2}$ bifurcations individually, we comment on $T_{1}$ and $T_{2}$ which involve mode interactions. The trajectory of IPO lies on $y=\Delta_{y}$ plane, hence bounded within a circular orbit of radius $\left(1-\Delta_{y}^{2}\right)^{1/2}-\alpha$ approximately, because the particle simply cannot penetrate the droplet. The SM mode hence suppresses the IPO mode due to the confinement; a larger $\Delta_{y}$ and/or $\alpha$ naturally shrinks the orbital displacement $\Delta_{xz}$ to be zero, when entering
$\mathbf{5}$  across $T_{1}$. On the contrary, the emergence of $T_{2}$ reflects the promotive effect of IPO on SM. Regarding this, we may surmise that when the particle starts orbiting on $y=0$ plane, it comes closer to the droplet interface; therefore, it suffers a greater compressional flow (as indicated by Fig. 5b) which results in stronger destabilizing viscous shear forces.

In summary, we have presented in this Letter, hydrodynamic-interaction-meditated dynamics of a particle inside a droplet in steady shear flow. We have numerically discovered several equilibrium solutions where the composite system exhibits spontaneous symmetry breaking and unsteady dynamics rising through supercritical pitchfork and Hopf bifurcations; the particle can execute spanwise migratory and/or in-plane orbital movement. The bifurcations are partially attributed to the droplet-induced disturbance flow characterized by a kinematically unstable stagnation point. We have performed a normal-form analysis to delineate the interplay between bifurcations, revealing the suppression of the Hopf bifurcation by migration and promotion of the pitchfork bifurcation by orbital motion. The interplay can be rationalized by the geometric confinement and the disturbance flow.

It is worth-pointing that the bifurcation dynamics might not be directly generalized to the two commonly adopted models of cells, capsule and vesicle featured with elastic membranes. The in-plane elastic stresses developed on the interface might considerably suppress the interior flow that influences the inclusion dynamics.

We envision that our results might potentially inspire new approaches of ’hydrodynamic centering’ composite systems like emulsions to obtain a uniform shell in addition to electric centering methods Bei et al. (2008); Tucker-Schwartz et al. (2010), or vice versa, using hydrodynamic effect to generate emulsions with pre-designed nonuniform shell thickness Hennequin et al. (2009) for programmed release of substances. We hope our study will motivate experiments in these directions. We plan to address in our future work the influences of non-uniform shear, geometric features and confinement of the setup, which are all relevant for practical applications.

The authors thank Jan Guzowski, Jérôme Hoepffner, Philippe Meliga, Arne Nordmark and Howard A. Stone for useful discussions. The computer time is provided by the Swiss National Supercomputing Centre (CSCS) under project ID s603 and by SNIC (Swedish National Infrastructure for Computing). A VR International Postdoc Grant from Swedish Research Council ’2015-06334’ (L.Z.) and an ERC starting grant ’SimCoMiCs 280117’ (F.G.) are gratefully acknowledged.