{'fraction_non_alphanumeric': 0.030855539971949508, 'fraction_numerical': 0.010813011808351808, 'mean_word_length': 4.699716348633316, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hydrodynamics and rheology of a vesicle doublet suspension

The dynamics of an adhesive vesicle doublet under various flow conditions is investigated numerically using a high-order, adaptive-in-time boundary integral method. In a quiescent flow, two nearby vesicles move slowly towards each other under the adhesive potential, pushing out fluid between them to form a vesicle doublet at equilibrium. A lubrication analysis on such draining of a thin film gives the dependencies of draining time on adhesion strength and separation distance that are in good agreement with numerical results. In a microfluid trap where the stagnation of an extensional flow is dynamically placed in the middle of a vesicle doublet through an active control loop, novel dynamics of a vesicle doublet is observed. Numerical simulations show that there exists a critical extensional flow rate above which adhesive interaction is overcome by the diverging stream, thus providing a simple method to measure the adhesion strength between two vesicle membranes. In a planar shear flow, numerical simulations reveal that a vesicle doublet may form provided that the adhesion strength is sufficiently large at a given vesicle reduced area. Once a doublet is formed, the effective shear viscosity of a dilute suspension of vesicle doublets is found to be a function of reduced area. Results from these numerical studies and analysis shed light on the hydrodynamic and rheological consequences of adhesive interactions between vesicles in a viscous fluid.


I. INTRODUCTION
Vesicles (closed fluid-filled phospholipid bilayer membranes) have been widely utilized as biological cell mimics in biophysics and material engineering [1][2][3]. The ever expanding applications of vesicles have encouraged detailed experimental, theoretical, and numerical investigations of vesicle dynamics in applied flow, electric fields, and magnetic fields. Experimental investigations of vesicles have uncovered vesicle properties and various novel applications [4][5][6]. Theoretical investigations are often limited to small-deformation analysis on a nearly spherical vesicle or spheroidal analysis on a spheroidal vesicle [7][8][9][10][11][12]. On the other hand, various numerical methods have been developed for simulating the transient hydrodynamics of vesicle suspensions [13][14][15][16].
Hydrodynamics of a single vesicle in Stokes flow has been extensively investigated. In a planar shear flow, the vesicle hydrodynamics is characterized by the reduced volume (reduced area in two dimensions), viscosity contrast between interior and exterior fluids, and shear rate of the imposed far-field fluid flow. In addition, a vesicle with a rigid particle inside is also investigated as a biological mimic of a eukaryotic cell with a nucleus that occupies nearly half of the intracellular volume [17]. Small-deformation analysis shows that a vesicle tank-treads in a planar shear flow for low viscosity contrast and shear rate. At high viscosity contrast the tank-treading dynamics transitions to tumbling dynamics [8,9], and this leads to a transition in the effective shear viscosity of the vesicle suspension [8,18] that is also validated by direct numerical simulations [19] and experiments [20][21][22]. Between tank-treading and tumbling vesicle dynamics, a breathing (tremble) mode is also observed [8,[23][24][25] to alter the effective shear viscosity. In an extensional flow (planar or uniaxial), vesicle shape dynamics depends sensitively on the vesicle reduced volume and the elastic capillary number [23,[26][27][28]: Asymmetric shape and oscillatory undulation of the vesicle membrane are two examples of the complex vesicle hydrodynamics in an extensional flow.
Collective hydrodynamics of vesicles is dictated by the vesicle-vesicle interactions. In a quiescent flow, vesicle-vesicle adhesion leads to the formation of vesicle doublets [29,30] or clusters [31,32]. As a model for red blood cell (RBC) aggregates, a simplified model for adhesive vesicle-vesicle interactions can reproduce vesicle shapes similar to those observed in experiments of fibrinogen-induced RBC aggregates [31,32]. Using the Lennard-Jones (L.-J.) potential for point-point interaction between two vesicle membranes without the molecular details of adhesive interactions between RBCs, Flormann et al. [32] found that, under strong adhesion, the membrane may buckle to a sigmoidal shape in the contact region. Such a sigmoidal vesicle shape is also observed in RBC doublets and explained in a slightly different model [30]. It remains unknown how adhesive interactions between vesicles might lead to different vesicle hydrodynamics that results in different rheological properties of a vesicle suspension. Such studies will provide useful insight into the rheology of blood flow and how to use nano solutes to control the hypertension by adjusting the adhesive interactions between RBCs.
In numerical simulations of a vesicle suspension, the adhesive interaction between two vesicles is often ignored [15,33] for numerical convenience, mostly because a small time step is often needed to resolve the dynamics of the lubrication thin film between two membranes with adhesion. However, adhesive interactions between RBCs lead to RBC clusters that are expected to alter the rheological properties of the RBC suspension [31,34]. Hydrodynamic simulations of adhesive vesicles in a quiescent flow have revealed physical insights to the equilibrium shapes of RBC aggregates in both experiments [32] and theory [30]. The main goal of this work is to investigate the hydrodynamics of adhesive vesicles in flow conditions common in microfluidics such as a planar shear flow and an extensional flow.
The electrostatic nature of lipid molecules (a hydrophobic tail and a hydrophilic head with an electric dipole) complicates the interactions between a lipid bilayer membrane and another bilayer membrane [35][36][37][38] or a solid (such as a glass substrate or a nano particle). Adhesion between a vesicle and a solid has been extensively studied with more focus on the static equilibrium shapes [39][40][41][42][43][44][45][46][47][48] than the transient adhesion process [49][50][51]. Adhesive interactions between lipid membranes are essential in many biomedical, biological, and biophysical processes. For example, vesicle adhesion is crucial to initiate membrane fusion and fission in endocytosis, exocytosis, and the transport of small vesicles through membrane surfaces. In the absence of an external electric field and ions in the solvents, it is reported that the adhesive interactions between two membranes can be well approximated by the L.-J.-type potential [32], that consists of a long-range attraction component and a short-range repulsion component [37], and the combination of the two gives rise to an equilibrium distance where the interaction potential is at its minimum [37].
The strength of membrane-membrane adhesion can be estimated through the membrane contact angle at the edge of contact zone [52,53]. Alternatively, Evans and Metcalfe [35] used the micropipette to measure the adhesive force between two vesicles bound by strong adhesion. They were able to measure the reduction in the free energy per unit area of membrane-membrane contact formation due to the van der Waals' attraction. For the adhesive interaction between two lipid bilayer membranes in buffer solutions, a typical range for the energy density is between 1 to 10 µJ/m 2 , similar to the adhesion energy density of 3 µJ/m 2 between two RBCs [32]. For adhesion interaction between a vesicle and a substrate [42], the boundary between weak and strong adhesion is around 1 µJ/m 2 : the vesicle-substrate interaction is "strong" when the adhesion energy density is larger than 1 µJ/m 2 and "weak" when the energy density is less than 1 µJ/m 2 . The ratio of total adhesion energy to the bending energy (bending rigidity of lipid bilayer membranes ∼ 10 −19 J) gives a measure of vesicle deformation in the presence of adhesion [52].
In this work we focus on regimes where such ratio is of order one, between the weak adhesion (adhesion energy/bending energy 1) and strong adhesion (adhesion energy/bending energy 1) regimes. In this regime the membrane deformation may increase the "contact area" between two vesicles and enhance the adhesion effects. The equilibration of vesicle membranes in the strong adhesion regime has been well-studied and documented (see [32,47,52] and references therein). However it is unclear how the adhesion couples to the vesicle hydrodynamics in this regime. This paper aims to address this question with quantitative characterization in terms of physical parameters.
In a quiescent environment, sub-micron size vesicles are found to form a doublet due to their van der Waals' attractive interactions [52]. Due to the strong van der Waals' adhesive force, the vesicles are far from spherical shape and the membrane is almost flat in the "contact" region. Gires et al. used small-deformation analysis to investigate the hydrodynamic interactions between two vesicles in a planar shear flow with a long separation distance [54]. They found that the vesicle interaction could be either attractive or repulsive depending on the organization of the two vesicles relative to the shear flow. To the best of authors' knowledge, the effects of close-range vesicle adhesion on their hydrodynamics under an external flow have not been studied and quantified. The goal of this work is to use state-of-the-art boundary integral simulations to numerically investigate the dynamics of two vesicles in both planar shear flow and extensional flow for a wide range of vesicle shapes and adhesive strength and distance.
Boundary integral equation (BIE) approaches are well-suited for solving the low Reynolds flow problems considered here as they lead to reduction in dimensionality and achieve high-order accuracy even for moving geometry problems. When the vesicles adhere, one major issue for BIE solvers is to resolve the vesicle-vesicle hydrodynamic forces which become nearly singular. We use an interpolation-based quadrature rule [55] to maintain high-order accuracy for all vesicle-vesicle separation distances. To overcome the numerical stiffness induced by the membrane bending forces and to control the global error, we employ a second-order spectral deferred implicit-explicit adaptive time stepping scheme [56].
This paper is organized as follows. In Section II we formulate our model for two-dimensional vesicle hydrodynamics with adhesive interactions between membranes of distinct vesicles. We simulate the adhesion process of two vesicles in a quiescent flow in Section III, where we also present a simple lubrication model to estimate how long it takes for two nearby vesicles to reach the separation distance set by the adhesion potential. In Section IV we numerically investigate the hydrodynamics of a vesicle doublet in a fluid trap where the stagnation point is actively controlled to be at the center of the vesicle doublet. From these results we propose a novel application of the microfluidic fluid trap to probe the adhesion strength between lipid bilayer membranes in solutions. In Section V we investigate how two vesicles may form a doublet as they move toward each other under a planar shear flow. We map out the parameter regions for bound/unbound vesicles under a planar shear flow, and we also investigate the effects of adhesive interactions on the rheological properties of a dilute suspension of vesicle doublets. Finally in Section VI we discuss the implications of our results and point out a few potential future directions.

II. GOVERNING EQUATIONS
We consider a suspension of locally inextensible vesicles in an unbounded two-dimensional viscous fluid. For simplicity, we assume the fluid viscosity both inside and outside the vesicles is the same; however, incorporating viscosity contrast is rather straightforward in our numerical algorithm. While we focus on a vesicle doublet suspension for the rest of the paper, here we provide the most general formulation for multiple vesicles interacting with each other through both hydrodynamic and adhesive forces. Individual vesicles are denoted as γ j , j = 1, . . . , M , and they are parameterized in arclength as x j (s, t). The union of all vesicles is denoted by γ. Given a background velocity u ∞ , the governing (dimensionless) equations are conservation of momentum velocity continuity x s · u s = 0, x ∈ γ, local inextensibility T n = ξ, x ∈ γ, stress balance on membranes (1) where u is the fluid velocity field, p is the pressure, and the scaled viscosity µ = 1 inside the vesicles. Outside the vesicle µ = µ e /µ i with µ i (µ e ) the viscosity of the interior (exterior) fluid. We set the viscosity ratio µ e /µ i = 1 for the rest of the paper. T = −pI + µ ∇u + ∇u T is the stress tensor, T is the jump across the membrane, and ξ is the traction that is the sum of the bending, stretching, and adhesion forces defined in equation (8) of Appendix A. The bending force, arising from the Helfrich energy model, is Bx = −κ b x ssss , where the subscript s denotes derivative with respect to arclength s. κ b is the bending rigidity modulus which we set to be 1 for all examples. The stretching force is T σ = (σx s ) s , where the tension, σ, acts as a Lagrange multiplier to satisfy the local extensibility constraint. The resistance to bending and stretching are standard assumptions on vesicles. In this work, we include an adhesive force, Ax, that we now describe.
The adhesion potential φ(z) in equation (2). Distances less than δ are repulsive and distances greater than δ are attractive. Center: Each point on a vesicle is repelled or attracted by all points on the other vesicles. Right: The resulting adhesive force is found by integrating over all points on all other vesicles (equation (3)). The color is the dot product of the adhesion force with the vector joining the vesicles. Therefore, the green and yellow regions are repelled by the other vesicle while the blue regions are attracted to the other vesicle.
The adhesion potential for a vesicle membrane has been modeled by a L.-J. type potential where H is the Hamaker constant and δ is the adhesion length scale, and z is the distance between a patch on the vesicle membrane and a point on the other object, which could be another vesicle membrane or a flat solid wall. The characteristics of the L.-J. potential are summarized in the left panel of Figure 1: The adhesive force is zero at the equilibrium distance z = δ. For large distances (z > δ) the adhesive potential is attractive while for small distances (z < δ) the potential is repulsive to prevent physical contact. The exponents (m, n) in equation (2) depend on the geometry and the molecular details of the two objects under adhesion [37]: (m, n) = (4, 2) corresponds to two flat, planar surfaces interacting with each other, and has been a common choice for studying membrane-solid adhesion [39-41, 50, 51, 57]. On the other hand, (m, n) = (12,6) corresponds to the L.-J. potential between two molecules, and has been used to model membrane-membrane adhesion [32]. In our case, the adhesion potential is between two small patches of lipid bilayer membrane (because the lipid molecules are coarse-grained in the continuum modeling). Thus a reasonable choice for (m, n) between two coarsegrained membrane patches would be between those for two planar surfaces and two point molecules.
Large value of m corresponds to a sharp increase in the repulsion force as two objects are within the separation distance δ. This poses a numerical challenge since the problem becomes stiff (i.e., requires a very small time step) for large m. The adaptive time-stepping BIE scheme makes it possible to simulate vesicle adhesive dynamics with specified numerical precision for reasonable computation time. We explore several combinations of (m, n) in the simulations of two vesicles forming a doublet in a quiescent flow and in a shear flow. We found that, as long as the close-range interaction is well resolved both in space (using the near-singular evaluations) and time (using the high-order adaptive time integration), there is very little difference in both the dynamic evolution and equilibrium configuration between (m, n) = (8, 6) and (m, n) = (4, 2). Therefore, in this work, we use (m, n) = (4, 2) to regulate the numerical stiffness introduced by the adhesive force.
Focusing on intermediate adhesive strengths, we use Hamaker constants ranging from O(10 −1 ) to O(10 0 ) times the bending rigidity modulus. We assume that the adhesion force from equation 2 applies between all pairs of points on different vesicles, and in Appendix B the net adhesive force at a point x on vesicle j is shown to be This summation of adhesive forces is illustrated in the middle panel of Figure 1. The right panel of Figure 1 is an example of the calculated adhesive force projected onto their center-of-mass vector. We notice that the adhesive force is repulsive in the contact region, while for the rest of the membranes there is an attractive force that acts to pull the two vesicles together. The summary of boundary integral formulation can be found in Appendix A.

III. ADHESION OF TWO VESICLES IN A QUIESCENT FLOW
We consider two identical vesicles suspended in a quiescent flow. Without any external forcing such as an imposed electric field or a fluid flow in the far-field, the long-range attraction pulls the vesicles together until their separation distance is close to δ. The L.-J. type potential prohibits physical contact between the two membranes and instead keeps them close to the separation distance δ. In Section III A, we compute the expected time for the vesicles to reach an equilibrium configuration and compare the analysis with numerical simulations. In Section III B, we characterize the effects of the adhesion parameters and the vesicles' reduced area on the shape of the contact region and the bending and adhesive energies.

A. Effects of the adhesion parameters on draining times
When two vesicles move towards (or away from) each other under a constant force F without any imposed external flow, the height, h, of the thin liquid film between two vesicles follows the draining dynamics [58] where K a is the area expansion modulus, µ is the viscosity of the exterior fluid, and R 0 is the radius of the undeformed vesicle. For the case of a constant forcing F (independent of separation distance h and time t), equation (4) can be easily integrated to relate the film thickness h and time t: where h(0) is the vesicle separation at t = 0. We note that F < (>)0 for an attractive (repulsive) interaction between two vesicles, and consequently h(t) decreases (increases) from the initial separation distance h(0). When the force is attractive, h decreases monotonically, and t in the above equation is the draining time that diverges as h → 0. When the force on each vesicle is a function of the separation distance h of the form with integers m > n ≥ 2. Attraction is dominant at "large" distances (h > δ) while repulsion is dominant at "small" distances (h < δ). Integration of equation (4) with the F (h) as defined in equation (5) gives the relationship between t and h. In this case, the relationship involves an integral of a function of the dimensionless variableh ≡ h/δ: and equation (6) tells us that both the adhesion strength H and the separation distance δ affect the time it takes for a pair of vesicles to reach the separation distance δ under adhesion force F in equation (5). In particular, the duration t is proportional to (i) the separation distance δ −5/3 (thus t → ∞ as δ → 0), and (ii) the adhesion strength H −1/3 (as in the constant forcing case [58]). From the above analysis we also expect that the scaling of t with respect to δ depends on the adhesion potential, while the scaling with respect to adhesion strength is independent of the specific form of the potential. To test the scaling of the duration t with respect to H and δ, we simulate the vesicle adhesion dynamics in a quiescent flow. Starting with two vesicles at a distance of twice the vesicle radius, the long-range attraction pulls the vesicles together. The left plot of Figure 2 shows the scaled adhesion energy versus time with H = 0.1 and three values of δ as labeled. Adhesion energy reaches minimum at equilibrium, and the scaled E adh evolves towards zero at equilibrium. The arrows indicate the times when the scaled adhesion energy reaches within 1% of equilibrium: t ∼ 96 for δ = 0.2, t ∼ 45 for δ = 0.3, and t ∼ 30 for δ = 0.4: The scaling with respect to adhesion strength H is illustrated in the right plot of Figure 2, where δ = 0.4 and H varies from 0.1 to 0.5 as labeled. Again the arrows indicate the times when the equilibrium is reached within 1%: From the above results, we conclude that the scaling in equation (6) captures the adhesion dynamics of two vesicles that interact with each other via the adhesion force in equation (5).

B. Effects of the adhesion parameters on equilibrium configuration of a vesicle pair
We again simulate two identical adhering vesicles with the adhesion force in equation (5). The initial vesicle separation is smaller than in the previous simulations so that the equilibrium configuration is achieved in a shorter time horizon. Once a vesicle doublet is formed, the membrane shape in the contact region depends on the adhesion strength relative to the membrane bending rigidity. For weak to moderate adhesion strength, vesicle membranes are flattened in the contact region while the rest of vesicle maintains a nearly spherical shape [35][36][37]52]. Under a strong adhesion, however, the vesicle membranes in the contact region buckle and form a sigmoidal shape that has also been observed in RBC doublets [29,30,32]. An external electric field is also able to buckle a vesicle membrane that is adhered to a solid substrate [47].
In this work we focus on weakly adhesive vesicles whose equilibrium shapes in a quiescent flow are shown in the left plot of Figure 3 for three reduced areas: ∆A = 0.95, 0.75, and 0.6 with H = 5 and δ = 0.2. The equilibrium vesicle shape for ∆A = 0.95 is the circular cap with a flat contact region. This is similar to the observed shapes of two vesicles under strong adhesive interaction in [52]. When vesicles are more deflated with a reduced area ∆A = 0.75, the equilibrium vesicle shape is elongated with a bigger contact region. This is consistent with the equilibrium shapes of a vesicle doublet under a (m, n) = (12, 6) L.-J. potential [32]. As the reduced area decreases further, we observe undulation of the vesicle membrane on the non-contact side while the contact region remains flat. The color coding along each curve is the tension of the vesicle membrane. We observe that the membrane tension in the contact region is very negative, indicating a dominant compression of membrane when the adhesion force is strong to keep the vesicles bound together.
The right plot of Figure 3 is a zoom (not to-scale) on the membranes in the contact region, where the membranes are not perfectly flat at all three values of reduced area: We observe that the membranes are slightly curved with a dip at the edge. Such membrane undulation in the contact region is predicted by lubrication analyses on an elastic membrane under adhesion with a solid substrate [51,57]. Results from the lubrication analysis show that this dip and slightly curved shape in the contact region are independent of the adhesion strength. The high membrane curvature at the edge may pose a problem for using the contact angle there to estimate the adhesion strength. This inspires us to investigate the possibility of using a dynamic fluid trap to measure the adhesion strength (see Section IV).
Also the dip at the edge of contact region is related to (but different from) the buckling of membrane under strong adhesion: A closer inspection on the dip at the edge shows that the membrane distance is smaller than the neutral separation distance δ there for more deflated vesicles. For A = 0.6 and 0.75 the adhesion force is mostly attractive in the contact region except at the edge, where the adhesive force turns repulsive. For A = 0.95, however, the adhesion force is mostly repulsive, leading to large tension in the rest of the vesicle membrane.
The left plot of Figure 4 shows the total bending energy and adhesion energy as a function of the reduced area. The smaller the reduced area, the more vesicle area is available for deformation and thus the bending energy is higher. In contrast, the total adhesion energy becomes more negative as the reduced area decreases. The sum of two energies is plotted in the right plot of Figure 4, where a local minimum in the total energy is found around ∆A = 0.85. Figure 5 demonstrates how the adhesion strength H and separation distance δ affect the equilibrium configuration of two vesicles under adhesive interactions. Both vesicles have a length of 2π and a reduced area of ∆A = 0.95. The total bending (left) and adhesion (right) energies at equilibrium are plotted against H for three values of the separation distance δ. We observe that for ∆A = 0.95 the equilibrium vesicle shape does not vary much with the adhesion strength H, while the total adhesion energy varies linearly with H.

IV. ADHESION OF TWO VESICLES IN AN EXTENSIONAL FLOW
The hydrodynamics of a single vesicle in an extensional flow has revealed novel nonlinear vesicle dynamics not found for a viscous drop [23,[26][27][28]. The stagnation point (zero-velocity at the origin) in an extensional flow, u = χ(−x, y), is a saddle point where two laminar streams converge along an axis and then diverge in the orthogonal direction. With an active control algorithm to place the stagnation point at a desirable location by adjusting the streaming flow strength with a feedback loop [59], a particle can be trapped at the stagnation point for long time scales to facilitate image acquisition or other detailed measures such as particle image velocimetry of flow inside and around the particle. To explore the application of such a fluid trap to measure the adhesion strength between two bound vesicle membranes, we propose the following experiment. Beginning with two identical vesicles in an equilibrium configuration that form a doublet with a flat contact region (Figure 3), we turn on the fluid trap with the stagnation point placed at the center of the vesicle doublet. Depending on the Hamaker constant, we expect that either the flow overcomes the adhesive force and the doublet is broken, or the adhesive force is sufficiently strong and the doublet reaches a stable stationary configuration.
Based on how the contact region aligns with the extensional flow, we can define the inclination angle of the vesicle doublet as illustrated in Figure 6. When θ = π, the doublet long axis is aligned with the far-field converging flow, and the diverging flow pulls the vesicles apart from each other and the long-range attractive force is essential to keep the two vesicles from separating. In contrast, when θ = π/2, the doublet long axis is aligned with the diverging flow, the converging flow pushes the vesicles towards each other and the short-range repulsive force is essential to keep the two vesicles at the separation distance. With the converging flow pushing the two vesicles towards the stagnation θ = π θ θ θ = π θ FIG. 6. The definition of the inclination angle of a vesicle doublet at the center of an extensional flow. An inclination angle of π indicates a doublet whose long axis is orthogonal to the diverging direction, while an inclination angle of π/2 indicates a doublet whose short axis is orthogonal to the diverging direction. The square mark is a stagnation point and the round point is located at the center of the vesicle. The initial configurations in this section all have an inclination angle of θ = π.
point dynamically placed at the doublet center, this configuration is expected to be more stable than the θ = π configuration in the left of Figure 6. This is because when θ = π the diverging flow is pulling the vesicles away from each other as shown in the illustration. However, when the doublet inclination angle is fixed at θ = π/2, the vesicle doublet stays intact if the adhesive force is strong enough to dominate the fluid drag. However, it is possible that the diverging flow may overcome adhesion and break up the doublet. It is also reasonable to expect a vesicle doublet to rotate from the less stable configuration θ = π to the more stable configuration θ = π/2, so that the fluid flow is pushing them together rather than pulling them apart. In the following numerical experiments we initially place a vesicle doublet with θ = π, so that the diverging flow may be strong enough to pull vesicles away from the doublet. At low extensional rates, we expect the doublet to stay bound at the fluid trap stagnation point. On the other hand, the vesicle doublet may become unstable and eventually separate at higher extensional rates. Thus we expect there to exist a critical extensional rate χ c above which the vesicle doublet cannot stay bound under a given adhesion potential. Therefore, the dependence of the critical extensional rate χ c on the adhesion potential and the mechanical vesicle properties provides a means to probe We observe that the doublet rotates from θ = π configuration towards the more stable θ = π/2 configuration, as shown in Figures 7. We also observe that with moderate extensional rates, the doublet falls short of aligning their long axis with the diverging direction. Moreover, the final inclination angle is closest to the stable θ = π/2 for the smallest extensional rates. While this doublet remains bound for all the considered extensional rates, a doublet with a reduced area of ∆A = 0.8 is split with a critical extensional rate χ c = 0.1 (Figure 8).
In Figure 9, we plot the inclination angle of the doublet as a function of time for four different extensional rates (left) and the final inclination angle for all doublets that reach an equilibrium state (middle). We observe that not only do smaller extensional rates result in smaller inclination angles, but smaller reduced areas also result in smaller inclination angles. We summarize the final inclination angle of the doublet in the right plot of Figure 9. The size of the round dots are scaled to the final inclination angle as defined in Figure 6. At smaller extensional rates, the vesicles come closer to aligning their long axis with the diverging direction. When the doublet is broken at reduced area ∆A = 0.80, the extensional rate is sufficiently large to align the long axis of the vesicle with the diverging direction and then the vesicles separate ( Figure 8). For very low extensional rates, the time horizon is insufficient for the doublet to reach an equilibrium state, and these simulations are marked with a square. The right plot of Figure 9 also summarizes the reduced areas and extensional rates that result in a bound vesicle doublet in a fluid trap: Parameter values with a blue mark result in a fluid trap and parameter values with a red mark result in vesicle separation. Finally, to further characterize the rotating dynamics of a doublet in the fluid trap, we compute and plot the vorticity and total force on a doublet with identical vesicles with reduced area ∆A = 0.7, Hamaker constant H = 0.7, separation distance δ = 0.4, and extensional rate χ = 0.7 in Figure 10. Starting from an initial inclination angle θ = π, the doublet starts to rotate around t ∼ 250. We notice that the net force on each vesicle is almost anti-parallel to the diverging flow in the far-field so to balance the hydrodynamic drag force that pulls the vesicles apart. Then, as the doublet starts to rotate, each vesicle experiences a drag force from the diverging streaming flows. For small extensional rate, the rotation continues until θ = π/2 so the doublet remains perfectly symmetric with respect to their mid-plane, and each vesicle experiences drag forces (from the diverging flow) that cancel each other. For slightly higher extensional rate in Figure 10, we observe that the two vesicles slide against each other (around t = 280) as the doublet rotates. The slight sliding motion leads to a transition in the contact region membrane shape from a bulging profile to an S-shape profile of nearly uniform height except near the edge of contact region. This adjustment greatly reduces the net force on each vesicle, leading to an equilibrium configuration with an inclination angle θ > π/2 as shown in the left of Figure 9. Finally, the vorticity and the flow between the vesicles approaches zero and the result is a fluid trap containing the doublet.
FIG. 10. The vorticity and total force of two vesicles in a doublet. The reduced area is ∆A = 0.7, the separation distance is δ = 0.4, the Hamaker constant is H = 0.7, and the extensional rate is χ = 0.7.

V. ADHESION OF TWO VESICLES IN A SHEAR FLOW
We consider two vesicles suspended in the shear flow u = χ(y, 0), where χ is the shear rate. The vesicles are placed on the y axis so that the background velocity drives the vesicles towards one another. In the absence of adhesion, once the vesicles are sufficiently close, they are deflected to opposite sides of the y axis, pass one another, and then separate. However, in the presence of adhesion, the vesicles can form a doublet for certain values of the separation distance δ, Hamaker constant H, shear rate χ, and reduced area ∆A. At all reduced areas, the vesicles have a length of 2π to be consistent with simulations in previous sections. Figure 11 shows snapshots of two vesicles that have formed a doublet and the color coding is the tension along the vesicle. The two vesicles are identical with reduced area ∆A = 0.9, shear rate χ = 0.5, separation distance δ = 0.4, and Hamaker constant H = 0.7. The vesicles in the doublet move in tandem rather than performing the tank-treading motion that is observed in the absence of adhesion. Similar to the quiescent example, the membrane tension is negative in the contact region indicating that the membrane is being compressed when the adhesive force is strongest.
We characterize the effect of the Hamaker constant and shear rate by determining a critical Hamaker constant, H c , which determines if the vesicles form a doublet or separate. We fix the separation distance to δ = 0.4 and the vesicle reduced area to ∆A = 0.9. In the left plot of Figure 12, we plot the minimum distance between the vesicles as a function of time for various Hamaker constants. The red curves correspond to Hamaker constants that are too small to bind the vesicles into a doublet, and the blue curves correspond to Hamaker constants that are sufficiently large to form a doublet.
Finally, we investigate the rheology of a suspension of a doublet by computing the effective viscosity of a doublet and compare it to the effective viscosity of a single tank-treading vesicle. The effective viscosity is defined as the viscosity of a homogeneous Newtonian fluid with the same energy dissipation per macroscopic element of fluid. In a simple shear flow, the intrinsic viscosity, [µ] is φ is the area fraction of vesicles, τ is the stress due to the vesicles, · is the spatial average, and ξ is the traction as defined in equation (8). In the left plot of Figure 13, we compare the intrinsic viscosity of a single tank-treading vesicle to a doublet. To validate our simulations, we superimpose (black marks) the intrinsic viscosity calculated by Ghigliotti et al. [19]. The presence of the doublet significantly increases the intrinsic viscosity at all the reduced areas. To further characterize the effect of adhesion, in the right plot of Figure 13 we decompose the intrinsic viscosity into the contributions from the bending and tension (blue), and the contribution from the adhesion (red). We also superimpose twice the intrinsic viscosity of a single tank-treading vesicle (dashed curve) to demonstrate that the bending and tension of the doublet behave similarly, but not identically, to a dilute suspensions of non-adhering tank-treading vesicles. We see that the effect of the adhesion on the intrinsic viscosity is largest for vesicles with small reduced areas. Left: The intrinsic viscosity of a tank-treading vesicle (blue) and a doublet (red). The shear rate is χ = 0.5, the Hamaker constant is H = 0.7, and the separation distance is δ = 0.4. The black marks denote intrinsic viscosity values computed by Ghigliotti et al. [19] (cf. Figure 5). Right: The decomposition of the intrinsic viscosity into the contributions for the bending and tension (blue) of the vesicles, and the contribution from the vesicle adhesion (red). Also included is twice the intrinsic viscosity of a single tank-treading vesicle (dashed black line).

VI. CONCLUSIONS
In this work we use a boundary integral formulation with adaptive time-stepping to simulate hydrodynamics of two vesicles with adhesive interactions. In a quiescent flow, two vesicles that are initially sufficiently far apart move towards each other under a long-range attraction. We use a lubrication theory to estimate the time required to reach the separation distance δ, and the theoretical scaling is in good agreement with numerical results. Once two membranes are within separation distance δ, the adhesive force turns repulsive and the membranes flatten to form a contact region. Our simulations show that the membranes in the contact region are actually curved with end points at the shortest distance. Once a vesicle doublet forms, we examine the dependence of membrane bending and adhesion energies on the reduced area, the Hamaker constant, and the separation distance.
Next we conduct a numerical experiment where a vesicle doublet is placed at the center of a fluid trap, which can be actively controlled in microfluidic channel so that fluid trap center is effectively the stagnation point of an extensional flow characterized by a extensional flow rate. At low flow rate, we find the doublet to rotate nearly ninety degrees to align with the flow such that the long axis of the doublet is parallel to the divergent axis and the convergent stream is pushing the two vesicles together. As the flow rate increases, the vesicle doublet rotates less, and when the flow rate exceeds the critical value, the diverging flow breaks the doublet structure by pulling the vesicles apart. These results indicate that it is possible to use the fluid trap to separate a vesicle doublet under adhesion, and thus provide a means to probe the adhesion strength between membranes. For a pair of µm-sized vesicles with a bending rigidity of ∼ 10 −19 J, an extensional flow rate of 0.5 s −1 is expected to separate a vesicle doublet with reduced area of ∆A = 0.8 and a Hamaker constant H = 0.7, which corresponds to ∼ 1 µJ/m 2 . These conditions are quite realizable in microfluidic experiments, and we hope that our simulations will inspire microfluidic experiments.
We also examine how adhesive interaction may dynamically lead to the formation of a vesicle doublet. We simulate two vesicles approaching each other in a planar shear flow, and examine how their adhesive interactions lead to doublet formation. Once a doublet forms, the two vesicle membranes rotate around each other as they deform dynamically. The usual tank-treading motion of a vesicle under shear flow is not observed in each of the two vesicles. We compute the effective shear viscosity of a dilute suspension of vesicle doublets, and found it to be more than twice the effective shear viscosity of a dilute suspension of single vesicles. Furthermore, we find that the membrane adhesion contribution to the shear viscosity increases with decreasing reduced area while the bending/tension contribution increases with the reduced area.
In our formulation we did not include any electrostatic interactions between membranes under adhesion. When the electrostatic interaction is important, an electro-osmotic pressure in the thin film is found to be responsible for the observed membrane undulation [47]. In addition, simulations presented in this work are for identical vesicles (same reduced area, length, and bending rigidity modulus) in the doublet. Flormann et al. demonstrated that asymmetric vesicle reduced area may lead to various equilibrium doublet shape such as male-female, asymmetric S-shape, and parachute shape. It is also possible that viscosity contrast may also lead to different equilibrium doublet shapes. Future work includes three dimensional simulations, dispersive vesicle properties (such as reduced area and bending modulus), viscosity contrast, and effects of confinement on adhesive interactions.
Flormann et al. studied the clustering and packing of an unbounded suspension of vesicles [32]. In engineering applications of vesicle emulsions, smaller vesicles are often enclosed in a big vesicle, as shown in Figure 14, where the left panel shows an initial condition of five vesicles suspended inside of a larger vesicle. The second panel from the left shows the membrane configuration and flow streamlines at an intermediate time. There is no imposed flow, so the dynamics is governed entirely by the vesicle forces, including adhesion. In the absence of confinement, we would expect the vesicles to align themselves as a rouleaux with the membranes forming a parachute, male-female, or S-shape configuration [32]. However, because of the bounding vesicle and the dense packing, the interior vesicles align themselves differently. The bounding membrane evolves dynamically together with all the interior vesicles to minimize the total membrane energy.
As a result of force balance, the equilibrium shape of the bounding vesicle is far from a symmetric shape while the interior smaller vesicles are of similar shape to each other. We also plot the tension (third) and adhesion (fourth) of the final equilibrium configuration. We remark that the adhesion force is attractive on all the vesicle membranes while large positive membrane tension is found only on the bounding vesicle. In particular, the interior vesicles form a flat contact region with a large negative tension. We are now conducting simulations of such vesicle configuration under a planar shear flow. Finally, it is not clear how thermal fluctuations may affect the hydrodynamics of vesicles under adhesion. For example, does the fluctuating hydrodynamics in the thin film between two vesicles enhance adhesion to keep vesicles bound under linear flows as speculated in cells [60]? Recently Liu et al. [61] used immersed boundary simulations to show that, at a separation distance of tens of nanometers, the thin film between the two membranes facilitate the coupling between membranes via strong hydrodynamic interactions. In particular, they demonstrate numerically that the fluctuation in one membrane is highly correlated to the other membrane without any physical contact. We are actively pursuing this direction with hydrodynamic modeling and simulations of adhesive membranes with thermal fluctuations. Perturbing x 1 tox 1 = x 1 + δx 1 , results in a new vesicleγ 1 , and the perturbed adhesive energy is φ( x 1 − x 2 )ds x2 dsx 1 , and the change in the energy is We now decompose the perturbation into normal and tangential components as δx 1 = y(s) = (ut + vn). The perturbed arclength term, to leading order, is where κ is the curvature. To leading order, inextensible perturbations satisfy u s + κv = 0, so the arclength term of γ 1 andγ 1 are identical to leading order. Therefore, ∇φ( x 1 − x 2 ) · yds x2 ds x1 and the adhesive force applied by vesicle 2 on vesicle 1 is γ2 ∇φ( x 1 − x 2 )ds x2 = −mHδ n γ k x − y x − y m+2 δ m−n − x − y m−n ds y .
When (m, n) = (4, 2), the above express becomes A similar expression holds for the adhesive force applied by vesicle 1 on vesicle 2, and equation (3) gives the adhesive force for a suspension of M vesicles.