{'fraction_non_alphanumeric': 0.042679787173179246, 'fraction_numerical': 0.014989415870473139, 'mean_word_length': 4.556261919898284, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bifurcation dynamics of a particle-encapsulating droplet in shear flow

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.

To understand the behavior of composite fluid particles such as nucleated cells and doubleemulsions 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.
Droplets, capsules and vesicles in flow often exhibit interestingly rich dynamics even in the linear shear flow [1][2][3][4][5][6][7][8][9]. 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 ∼ 80% of themselves in volume [10]; double-emulsions playing an important role in chemical and pharmaceutical engineering are featured with a core-shell geometry [11][12][13] ; 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 [14][15][16].
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 [17][18][19][20] all assumed their compound structures to be concentric, preserving the rotational symmetry of order 2 (C2) about the y axis and reflection symmetry about the y = 0 shear plane (see Fig. 1a). The symmetries do hold for a single sheardriven deformable particle which attains a steady ellipsoidal shape undergoing tank-treading motion [21][22][23]. 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 3D hydrodynamic simulations of a compound particle-droplet subjected to unbounded shear U ∞ = G · x, in the creep-ing flow regime, where the only non-zero component G xz =γ represents the shear rate (Fig. 1a). The incompressible Stokes equations are solved by a boundary integral method (see Supplemental material and [24] for details). The immiscible Newtonian fluids inside and outside the droplet have the same viscosity η; its surfactantfree interface has a uniform surface tension σ. 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 [25,26]. The radii of the particle and the undeformed droplet are a and R respectively; the size ratio is denoted by α = a/R with α ∈ (0, 1). The capillary number Ca = ηγR/σ 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 d = (d x , d y , d z ), then focusing on the time evolution of d y and d xz = d 2 representing the spanwise and in-plane displacements re-spectively. ∆ y and ∆ xz denote their equilibrium values when the system reaches a steady or time-periodic state. Evolution of the in-plane and spanwise displacements dxz (green) and dy (red) between the centers of particle (blue) and droplet (pink), for size ratio α = 0.4 and Ca = 0.2 (first row), 0.275 (second row) and 0.29 (third row). The second (resp. third) column shows the profile of the composite system at equilibrium, cut by x = 0 (resp. y = 0) plane. The particle position is shown at different instants within a period for Ca = 0.275, the insets of (d) and (f) display the limit cycle solution.
We show the evolution of displacements in Fig. 2, presenting three typical Ca-dependent scenarios for a particle of size ratio α = 0.4 (see Supplemental videos). When 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 ∆ y ≈ 0.48 indicating the broken reflection symmetry. The particle rotates steadily near the lateral edge of the droplet interface (Fig. 2b). Increasing Ca to 0.275, ∆ y decreases to 0.21 approximately, while d xz reaches a time-periodic equilibrium cycle with a maximum of ∆ max xz ≈ 0.11 and a minimum of ∆ min xz ≈ 0.04. The particle follows an orbital trajectory on the y = ∆ y plane as it reaches a limit cycle solution in the (x, z) space (Fig. 2d, f), implying that the C2 symmetry and time invariance are also broken. At Ca = 0.29, the system recovers steadiness and concentricity, ∆ y = ∆ xz = 0. These scenarios suggest the appearance of bifurcating solutions by reducing Ca: above a critical value Ca c (α), the composite system stays concentric, corresponding to a stable fixed point solution; it bifurcates across Ca c (α) towards a steady spanwise migration (SM) and/or in-plane orbiting (IPO) motion.
An investigation spanning the (Ca, α) space further re-veals 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∆ y and∆ min/max xz are shown in Fig. 3. For all α, we observe∆ y decreases with Ca, becoming zero when Ca exceeds a critical value Ca c (α), so does ∆ min/max xz . They both vary quadratically in the vicinity of their corresponding Ca c (α). This is confirmed by the linear fitting of∆ y and,Ã xz =∆ max xz −∆ min xz indicating the oscillating amplitude, versus [Ca c (α) − Ca] 1/2 (Fig. 3d), where Ca c (α) 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 Ca and α (Fig. 3). For α = 0.3, the in-plane amplitudes (cyan) reach their maxima around Ca = 0.3, from where they decrease almost linearly/quadratically with decreasing/increasing 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 α = 0.4, the spanwise offsets coincide precisely with the decoupled counterparts for all Ca, while the in-plane amplitudes exhibit non-monotonic Ca-dependence as for α = 0.3 and they are below the decoupled values. For α = 0.5, perfect coincidence between the spanwise offsets also holds as in the α = 0.4 case, while ∆ min/max xz ≡ 0 for all Ca, i.e., the Hopf bifurcation is inhibited.
The complexity is better unraveled by the parametric portrait quartering the (Ca, α) 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 (Ca * , α * ) ≈ (0.286, 0.4) is pinpointed at the intersection of the two marginal curves H ± 1 (circle) and P ± 2 (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 (Ca * , α * ) 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 bifurca-  tion in Ref. [27] where the amplitudes are independent of phase evolution. Denoting the square of the in-plane and spanwise amplitudes by ρ 1 = ∆ max xz − ∆ min xz 2 and ρ 2 = ∆ 2 y , the truncated amplitude system is expressed aṡ where µ 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 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 ρ 1 ρ 2 2 and ρ 2 ρ 2 1 without changing the bifurcation topology [27]. By introducing new phase variables ξ 1 = −p 11 ρ 1 and ξ 2 = −p 22 ρ 2 , we obtainξ where θ = p 12 /p 22 and δ = p 21 /p 11 . Applying at leading order the affine transformation in the vicinity of (Ca * , α * ), we map the parameter space from (Ca, α) to (µ 1 , µ 2 ) (inset of Fig. 4b), where k 1 ≈ −2.3 and k 2 ≈ −40 denote the slope of H 1 and P 2 curves at (Ca * , α * ), with C ≈ 0.043 derived from the growth rate of ρ i . The slopes of T 1 and T 2 further determine θ ≈ 0.71 and δ ≈ −0.69. The parametric portrait therefore corresponds to case III described in Ref. [27], characterized by six regions: 1 corresponding to 'concentric'; 2 to 'pure IPO'; 3 and 4 to 'mixed' separated by P + 2 ; 5 and 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 2 ∼ 5 but it is suppressed due to the nonlinear interaction with SM, as reflected by the phase portraits of 3 ∼ 5 all including an unstable saddle-node equilibrium (µ 1 , 0), as well as by the sign of θ. Consequently, pure IPO only survives in 2 . Besides, without IPO, pure SM is stable in regions 4 ∼ 6 , while IPO promotes SM to expand its locus further to 3 that indeed involves a stable equilibrium (µ 1 , µ 2 ). This promotion results from the sign of δ. 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 Ca=0.325 Ca=0.2 with imposed concentricity in U ∞ E (x), Fig. 5 displays the droplet-induced disturbance flowsû 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 = ±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 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 Ca = 0.1; while the former dominates the latter at Ca = 0.2, hence overcoming the kinematic instability and leading to a concentric preference. This might explain the quenching of IPO when 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 ±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 Ca > 0. Moreover, we conduct simulations fixing a certain spanwise offset d fix y > 0 with d x = d z = 0 for the shear flow, recording the spanwise hydrodynamic forces F Y = F p Y + F η Y on the particle (Fig. 5c), where F p Y (resp. F η Y ) represents the pressure (resp. viscous) contribution which centralizes (resp. decentralizes) the particle. As shown, the viscous shear force F η Y | s accounts for the major contribution to F η Y , supporting the above arguments of compressioninduced viscous destabilization. The pressure force becomes stronger with Ca and dominates the viscous part when Ca exceeds a critical value. In fact, the droplet with larger 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 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 = ∆ y plane, hence bounded within a circular orbit of radius 1 − ∆ 2 y 1/2 − α approximately, because the particle simply cannot penetrate the droplet. The SM mode hence suppresses the IPO mode due to the confinement; a larger ∆ y and/or α naturally shrinks the orbital displacement ∆ xz to be zero, when entering 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 [28,29], or vice versa, using hydrodynamic effect to generate emulsions with pre-designed nonuniform shell thickness [30] 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.