Dynamics of an elastoviscoplastic droplet in a Newtonian medium under shear flow

The dynamics of a single elastoviscoplastic (EVP) drop immersed in plane shear flow of a Newtonian fluid is studied by 3D direct numerical simulations using a finite-difference and level-set method combined with the Saramito model for the EVP fluid. This model gives rise to a yield stress behavior, where the unyielded state of the material is described as a Kelvin-Voigt viscoelastic solid and the yielded state as a viscoelastic Oldroyd-B fluid. Yielding of an initially solid drop of Carbopol is simulated under successively increasing shear rates. We proceed to examine the roles of nondimensional parameters on the yielding process; in particular, the Bingham number ($Bi$), capillary number ($Ca$), Weissenberg number ($Wi$), and ratio of solvent and total drop viscosity. We find that all of these parameters, and not only $Bi$, have a significant influence on the drop dynamics. Numerical simulations predict that the volume of the unyielded region inside the droplet increases with $Bi$ and $Wi$, while it decreases with $Ca$ at low $Wi$ and $Bi$. A new regime map is obtained for the prediction of the yielded, unyielded, and partly yielded modes as a function of $Bi$ and $Wi$. The drop deformation is studied and explained by examining the stresses in the vicinity of the drop interface. The deformation has a complex dependence on $Bi$ and $Wi$. At low $Bi$, the droplet deformation shows a nonmonotonic behavior with an increasing $Wi$. In contrast, at moderate and high $Bi$, droplet deformation always increases with $Wi$. Moreover, it is found that the deformation increases with $Ca$ and with the solvent to total drop viscosity ratio. A simple ordinary differential equation model is developed to explain the various behaviours observed numerically. The presented results are in contrast with the heuristic idea that viscoelasticity in the dispersed phase always inhibits deformation.


I. INTRODUCTION
This manuscript considers the time-dependent, three-dimensional dynamics of an elastoviscoplastic (non-Newtonian) droplet surrounded by a Newtonian fluid in a simple shear flow. The dynamics of droplets in shear flows have attracted significant interest, mainly because this flow case governs the behaviour of dilute emulsions. Emulsions in turn play a pivotal role in a variety of applications, for example in food industry, chemical processing and biological materials. In order to control these industrial processes, it is hence important to understand the dynamics of a droplet in shear flow 1 . However, most studies so far are restricted to the case where both the droplet and the other fluid are Newtonian. In real life, both components are often non-Newtonian. It has been shown that non-Newtonian fluids a) Electronic mail: outi@mech.kth.se exhibit exotic behaviors which significantly affect the drop deformation and breakup. This study considers a droplet of a particular type of non-Newtonian fluid, which exhibits yieldstress (plasticity) along with elasticity. Droplets of such fluids have not been numerically simulated before, to the authors' knowledge. In simple shear flow, no results for droplets of any kind of yield-stress fluids (including purely viscoplastic) could be found in the literature; our results at low Weissenberg numbers may shed light on the drop behaviour in the viscoplastic limit.
Yield-stress fluids behave like solids below a local stress threshold (the yield stress), and flow like liquids above this threshold 2 . Yield-stress fluids can be found in geophysical applications, such as mudslides and the tectonic dynamic of the Earth. They are also found in industrial applications such as mining operations, the conversion of biomass into fuel, and the petroleum industry, to name a few. Biological and smart materials can have a yield stress, making yield-stress fluid flows relevant for problems in physiology, biolocomotion, tissue engineering, and beyond. Most of these applications deal with multiphase flows or interfacial flows (Maleki and Hormozi 3 , Chaparian and Frigaard 4 , Maleki et al. 5 among many others). Therefore, there is a compelling need to study droplet-laden flows with yield stress and predict their flow dynamics in various situations, including three-dimensional and inertial flows.
The simplest physical model of a yield-stress fluid assumes that the material is a rigid solid at stresses lower than the yield stress. This assumption leads to purely viscoplastic models, such as the Bingham model 6 , where the yielded fluid flow is Newtonian, and the Herschel-Bulkley model 7 , where the yielded fluid flow is shear-thinning. However, many yield-stress fluids, of which Carbopol is one example, deform like elastic solids in the unyielded state and behave as viscoelastic liquids in the yielded state, displaying elastic (E), viscous (V) and plastic (P) properties. Numerical simulations of EVP fluid flows are not a straightforward task due to the inherent nonlinearity of the governing equations. Nevertheless, numerical simulations are needed to provide quantitative information that is extremely difficult to access experimentally in yield-stress fluids (for example, the yielded and unyielded regions, and velocity fields and stress fields separated into different contributions), and also understanding the physics of the interaction between droplets and the surrounding fluid.
Numerical simulations have already helped to reveal elastic effects in yield-stress fluid flows around particles, in liquid foams and Carbopol. First, Dollet and Graner 8 performed experimental measurements for the flow of liquid foam around a circular obstacle, where they observed an overshoot of the velocity (the so-called negative wake) behind the obstacle. Then, Cheddadi et al. 9 simulated the flow of an EVP fluid around a circular obstacle by employing Saramito's EVP model 10 . The numerical simulation using the EVP model captured the negative wake. A purely viscoplastic flow model (Bingham model) on the other hand always predicted fore-aft symmetry and the lack of a negative wake, in contrast with the aforementioned experimental observations. The numerical simulations could hence prove that the negative wake was the consequence of the foam elasticity. Recently, the loss of the fore-aft symmetry and the formation of the negative wake around a single particle sedimenting in a Carbopol solution was captured by transient numerical calculations by Fraggedakis et al. 11 by adopting the Saramito EVP model, in agreement with the experimental observations in Carbopol gel 12 . The elastic effects on viscoplastic fluid flows have also been addressed in numerical simulations of the EVP fluids through an axisymmetric expansion-contraction geometry 13 . It was observed that elasticity alters the shape and the position of the yield surface remarkably, and elasticity needs to be included to reach qualitative agreement with experimental observations for the flow of Carbopol aqueous solutions. More recently, elastoviscoplastic effects have been analysed in porous media flow 14 and in turbulent channel flow 15 .
Regarding droplets in non-Newtonian fluids, a large body of literature has been devoted to viscoelastic fluids. Here, we only review studies on droplets in low Reynolds number viscoelastic shear flows for conciseness, although it needs to be mentioned that elastic effects (loss of fore-aft symmetry at low Reynolds) are commonly seen for buoyancy-driven bubbles and droplets 16,17 . In simple shear flows, viscoelasticity in the drop tends to reduce drop de-formation. The steady state deformation of an Oldroyd-B viscoelastic droplet under shear in Newtonian fluid was considered in several works [18][19][20] . With increasing elasticity (increasing Weissenberg number), the droplet deformed less, while the inclination angle of the droplet increased. Similar behaviour was found by Mukherjee and Sarkar 21 in viscosity-matched systems. However, when the droplet viscosity was higher than that of the surrounding fluid, they found that the droplet deformation started to increase with elasticity, at high Weissenberg numbers. The stresses on a single droplet were analysed and generalized towards stresses in dilute emulsions by Aggarwal and Sarkar 22 . It was found that the effective shear viscosity of the emulsion was not much affected by droplet phase viscoelasticity, while the normal stress differences increased with drop viscoelasticity at high capillary numbers. The effect of confinement of the viscoelastic droplet in microfluidic applications has also been analysed 23,24 , whereby it was found that confinement increased the droplet deformation and promoted breakup, and the viscoelastic stresses also increased.
While droplet viscoelasticity often just reduces deformation at the steady state, the timedependent dynamics may change qualitatively. A new breakup mode due to droplet viscoelasticity was found by Verhulst et al. 25 . The effect of maximal extensibility of the polymers (using the FENE-P numerical model) was studied by Gupta and Sbragaglia 26 , who found a nontrivial relation between the polymer extensibility and the droplet breakup behaviour. In systems where the droplet viscosity was much higher than that of the surrounding fluid, the droplet was seen to oscillate periodically 21 .
Problems where the surrounding fluid (matrix) is viscoelastic but the droplet is Newtonian have also been analysed. As the present work deals with the opposite configuration, we will not attempt to review these results comprehensively. However, Verhulst et al. 27 compared the effects of droplet viscoelasticity and matrix viscoelasticity experimentally and computationally, and found that the matrix viscoelasticity generally has a much larger effect on the droplet deformation and orientation. Early works on the effect of Weissenberg (or Deborah) number gave contradictory results: in some experiments, the droplet deformation decreased with matrix elasticity 28 , while other experiments found that the deformation increased 29,30 . The main reason for this apparent contradiction is that matrix viscoelasticity has a nonmonotoneous effect when the elasticity is increased: it reduces the droplet deformation at low Weissenberg numbers (weak viscoelasticity), while it increases the deformation at high Weissenberg numbers (strong viscoelasticity). The non-monotoneous dependence on the Weissenberg number was first discovered in the computational study of Yue et al. 18 , and later in a fully three-dimensional simulation by Aggarwal and Sarkar 20 .
The present work investigates how the dynamics of the droplet in shear flow change, when the droplet is elastoviscoplastic rather than viscoelastic. The EVP behavior is obtained from the Saramito 10 constitutive equation, without shear-thinning or time-dependent rheology effects. The manuscript is organized as follows. The governing equations and their numerical solution are described in Section II. In Section III, the flow case and its characteristic nondimensional parameters are presented. Numerical results are shown in Section IV for varying Bingham, Weissenberg and capillary numbers, and solvent-to-total viscosity ratio. We consider the distribution of yielded and unyielded regions inside the droplet. Next in Section V, we look at how the droplet deformation and inclination angle change when the above parameters are varied, and explain the change from a simple ordinary differential equation (ODE) model for fully unyielded droplets. Finally, we study the stress distributions along the interface (Section V A 1). The work is concluded in Section VI. Finally, the formulation of the ODE model is included as an Appendix A.

II. GOVERNING EQUATIONS AND THEIR NUMERICAL SOLUTION
The governing equations for this work are the momentum balance equations for two incompressible, immiscible fluids, combined with the constitutive equation for the elastoviscoplastic stresses. The numerical method that allows 3D numerical simulations of these equations is described in Izbassarov et al. 31 . At the interface between the two fluids, veloc-ity and tangential stress shall remain continuous, while there is a jump in the normal stress due to surface tension. These conditions were explicitly formulated in Tammisola et al. 32 , where a two-domain approach was used. In the present work however, the coupling terms between the two fluids (droplet and matrix) are treated using the level set method 33 . Hence, we adopt a continuous formulation of the governing equations 31,34 , with a delta function representing the normal stress jump due to surface tension. The level-set method used here is similar to the one used in our other works with droplets in low Reynolds number flows 35,36 , with the difference that the ghost-fluid method (GFM) is here replaced by a more diffusive continuum surface force (CSF) method, to simplify the treatment of the non-Newtonian stresses.
In the following, all dimensional quantities are denoted by stars. The governing equations for the dynamics of an incompressible flow of two immiscible fluids can be written as: where u * = u * (x * , t * ) is the velocity field, p * = p * (x * , t * ) the pressure field and τ * = τ * (x * , t * ) an extra stress tensor defined below. The function φ * (x * , t * ) is the level-set function approximating the signed distance from the interface. Hence, φ * = 0 denotes the interface, φ * > 0 denotes fluid 1 and φ * < 0 fluid 2. The last term in Eq. (1b) is a body force due to surface tension, where n * is the unit normal vector to the interface, κ * the local mean curvature, δ a regularized delta function and Γ * the surface tension. Finally, ρ * and µ * s are the density and solvent viscosity of the fluid. In the present study, the elastoviscoplastic (EVP) effects in the flow are reproduced by the extra stress tensor τ * , defined with the elastoviscoplastic Saramito model as where λ * is the relaxation time, µ * p is polymeric viscosity, τ * 0 is the yield stress, τ * d is the deviatoric part of τ * and | · | denotes the usual Euclidean matrix norm. The second deviatoric stress tensor τ * d and its magnitude is defined as in previous works 11,14,15,31 : We can observe that for τ * 0 = 0, the Oldroyd-B polymeric liquid model is recovered, and that if yield stress is large (τ * 0 → ∞), the equation describes a viscoelastic solid. The density, the solvent and the polymeric viscosities and the relaxation time vary across the fluid interface and are expressed in a continuous form as where φ is the level set function denoting the distance to the interface, the subscripts 1 and 2 denote the properties of the bulk and drop fluids, respectively, and H(φ) is the Heaviside function.
The flow equations (1a & 1b) are solved in parallel, fully coupled with the EVP model equations (2), by the numerical method recently presented in Izbassarov et al. 31 . The flow and the EVP model equations are solved on a staggered uniform Cartesian grid using a projection method. The spatial derivatives are approximated using second-order central differences, apart from the advection terms in the EVP constitutive, level-set advection and reinitialization equations, which are discretized by the fifth-order weighted essentially non-oscillatory (WENO) scheme 37 . The time integration is performed by the second-order Adams-Bashforth scheme (AB2) for both flow and elastoviscoplastic model equations. The time-integration in the level set advection and reinitialization equations is performed using a three-stage total-variation-diminishing (TVD) third-order Runge-Kutta scheme 38 . A complete description of the level-set method can be found in Sussman et al. 33 and Ge et al. 35 , and the treatment of the EVP multiphase flow in Izbassarov et al. 31 .

III. PROBLEM STATEMENT
In the present work, we consider a three-dimensional elastoviscoplastic (EVP) droplet in a simple shear flow of a Newtonian liquid. The configuration is sketched in figure 1, and x * is the streamwise, y * the vertical and z * the spanwise coordinate. The droplet is initially spherical and fully solid with a radius R * . Velocities U * and −U * in the streamwise x * −direction are imposed at the upper and lower walls, respectively, to obtain a shear ratė γ * = 2U * /L * y , which deforms the droplet. Under shear flow, the droplet may yield partly or fully, i.e. go from solid to liquid. The solid region is represented by the black colour in the following. As boundary conditions, we impose no-slip at the walls, and periodic boundary conditions in the streamwise and spanwise directions. The channel height is L * y = 8R * , resulting in a weak confinement of the droplet.
The lengths are nondimensionalized with R * , and velocities withγ * R * (whereby the shear flow time scaleγ * −1 has been chosen). The stresses and the pressure are nondimensionalized by the viscous stress scale of the outer flow µ * 1γ * . The outer fluid (i.e. matrix) density ρ * 1 and dynamic viscosity µ * 1 are taken as the reference values for nondimensionalisation. The problem is characterized by six non-dimensional parameters. Firstly, we have the droplet Reynolds number Re = ρ * 1γ * R * 2 /µ * 1 , and the capillary number Ca = R * γ * µ * 1 /Γ * characterizing the importance of surface tension compared to viscous forces. Furthermore, because the droplet is elastoviscoplastic, the Bingham number Bi = τ * 0 /(µ * 1γ * ) emerges as a ratio between the droplet yield stress τ * 0 and the matrix viscous stress. The elastic effects inside the droplet are characterized by the Weissenberg number W i = λ * γ * , where λ * is the elastic time scale. The last two parameters are the droplet-to-matrix viscosity ratio k µ = µ * 2 /µ * 1 and the density ratio k ρ = ρ * 2 /ρ * 1 . In the previous definitions, µ * 2 = µ * s,2 + µ * p,2 and µ * s,2 are the total (sum of the solvent and polymeric viscosities) and the solvent viscosity of the fluid 2, respectively. The solvent viscosity ratio is defined as β = µ * s,2 /µ * 2 . Considering that only the droplet (fluid 2) is EVP, the equations 1 & 2 can be rewritten in a nondimensional form as where ρ and µ s are the nondimensional counterparts of ρ * and µ * s in the previous section. To quantify the deformation of the droplet in the shear plane (x − y plane), we use the Taylor deformation parameter D = (L − B)/(L + B), where L and B are the major and minor axis of the equivalent ellipsoid in the middle plane. The inclination angle of the drop axis with respect to the streamwise direction is denoted by θ. Using a methodology developed by Ramanujan and Pozrikidis 39 , the equivalent ellipsoid is obtained by the inertia tensor of the drop.
In the following, numerical simulations are performed to study the dynamics of EVP droplet in a Newtonian fluid under shear flow. The computational domain (see figure 1) has dimensions L x ×L y ×L z = 16R×8R×8R and is resolved by a grid of ∆x = ∆y = ∆z = R/16 (32 points per drop diameter) in all the results presented in this paper. A grid convergence study has been performed ( Fig. 2) to ensure that the solutions are grid independent, i.e., the spatial error is below 3%. The solver has also been validated against the results of Verhulst et al. 27 for a viscoelastic droplet in a Newtonian shear flow (Fig. 2), and several other validation viscoelastic and EVP validation cases can be found in Izbassarov et al. 31 .   (Table I). The unyielded region is obtained from

IV. RESULTS ON THE EVP DROPLET YIELDING PROCESS
A. Droplet of a yield-stress fluid under shear -a numerical experiment We start the study by imitating an experiment, where a droplet of a given EVP material is placed in simple shear flow, under different levels of shear (γ * = U * /H * ). The EVP material properties are based on the data provided by 11 and are given in Table I & II. It is same Carbopol solution, but all of its parameters vary with shear rate. Note that here we employed a low but numerically tractable value of β, while the actual value for Carbopol is lower. Also for numerical reasons, we increased the minimal Reynolds number to Re = 0.05. Figure 3 shows the droplet interface at steady state, together with the unyielded regions defined by F = 10 −3 , where F = 1 − (Bi/|τ d |), at differentγ * = 0.32, 0.43 and 0.64s −1 . We have observed that the shape of the unyielded region is insensitive to smaller values of F, and hence it is fixed for all results presented in this work. At the lowest shear rate, the droplet remains solid almost everywhere. Similarly to an experiment, while keeping the material constant, we now increase the shear rate. Increasing the shear rate changes its state from fully unyielded to partly yielded and finally to fully yielded (fluidized everywhere). This qualitative development is expected for any EVP droplet because its Bingham number decreases with shear; the shear stress imposed by the external shear flow becomes stronger than the droplet yield stress.
In this dimensional experiment, W i, Bi and Ca all change simultaneously, and therefore, the effects of each parameter on the yielding process are difficult to separate and understand. In the next sections, we will adopt a nondimensional formulation and study the effects of each nondimensional parameter in detail.

B. Time evolution of the yielding process
We will start by characterizing the flow in and around an EVP droplet in time. Following Ref. 27 , 40 , the base case parameters are chosen as Re = 0.05, Ca = 0.32, W i = 2.31, β = 0.68, k ρ = 1 and k µ = 1.5. The parameters are fixed to the base case values unless indicated otherwise. The droplet has initially a spherical shape, and we follow its time evolution under suddenly imposed constant shear (constant plate velocity). This process differs considerably depending on the droplet elasticity, which will be demonstrated by comparing the process two different Weissenberg numbers: W i = 0.01 (nearly viscoplastic droplet) and W i = 0.5 (finite elasticity).
The time-evolution of the solid region for fully unyielded case at Bi = 0.8 and W i = 0.01 is shown in Fig. 4. As can be seen in the figure, the near-viscoplastic droplet yields fully and fast (by t = 0.0189), first remaining spherical. Eventually at t > 3, the droplet has deformed by the shear flow and attained a nearly ellipsoidal shape. This result is intuitive, because the near-viscoplastic droplet behaves like a rigid solid and cannot deform unless fluidized. This leads to the observed steady states for near-viscoplastic droplets throughout the paper: they either contain solid regions at steady state but have not deformed from the spherical shape, or have fluidized completely and deformed to an elliptical shape. This can be compared to the yielding process of a droplet at the same Bingham number Bi = 0.8 but at a larger Weissenberg number W i = 0.5, as depicted in Fig 5. In contrast to the near-viscoplastic case, the unyielding process is slow. The droplet first yields in the vicinity of the interface (t = 3.15), then solid region starts to oscillate yielding and unyielding diagonally across the droplet (3.15< t <54.6), until it finally attains a steady state with two separate solid regions (t = 54.6). The material yields fast at low W i, whereas at larger W i the material is more elastic, leading to a slower yielding process and a larger solid region. Elasticity enables deformation at the solid state, which allows elastic and capillary stresses to redistribute inside the droplet and counteract the shear force in some regions.
Solid (unyielded) regions in yield-stress fluids can either be stagnant (with zero velocity), such as fouling zones in porous media 41 , or they can move as solid plugs. In order to better understand the flow physics inside the droplet, where solid and fluidized regions may coexist, we show the velocity vectors at equilibrium for the more elastic EVP droplet shown before (Bi = 0.8 and W i = 0.5), along with a near-viscoplastic droplet in its completely unyielded state (Bi = 8 and W i = 0.01) in Fig. 6, right column. Clearly, in both cases the droplet experiences a tank-treading motion similarly to capsules and elastic particles in simple shear flow. By comparing the velocity distributions with the unyielded zones shown in the left column, we observe moving unyielded zones for both regimes, i.e. the material keeps moving inside the droplet. For the near-viscoplastic droplet, the movement is a simple solid body rotation. For the more elastic EVP droplet, the material continuosly unyields and yields when passing through boundaries of solid zones, while the solid zones remain at same place in steady state. This shows that the solid regions of the droplet are pseudoplugs, which have been observed in EVP flow in porous media 14,41 . We can also observe that the velocity vectors do not show any change in direction or character when entering or leaving the pseudoplug. This may seem unexpected, but is typical for pseudoplugs in yield-stress fluids. Pseudoplugs in cavity flows and cylinder flows of a Bingham fluid are shown in 42 , together with streamlines of the velocity field. In all those examples, the velocity field seems largely unaffected by the presence of pseudoplugs.
It is worth noting that the tank-treading motion constrains the maximum deformation of a material element inside the droplet. This is fundamentally different from parallel shear flows, where infinite deformation is possible. The deformation constraint may partly account for the surprising yielding behaviours at increasing Weissenberg numbers, discussed in the next section.

C. Effect of Weissenberg and capillary numbers
The dynamics of an elastoviscoplastic drop is governed by a subtle interplay between various competing effects, and we move further to explore their respective roles. First, we look at how the distribution of fluid (yielded) and solid (unyielded) regions inside the droplet changes with the Bingham number and the Weissenberg number, all other parameters kept constant. Figure 7 shows the droplet interface at steady state, together with the unyielded regions, at twelve different parameter combinations of Bi and W i. It is worth remembering that the left column contains only near-viscoplastic droplets and that elasticity increases to the right in the figure.
The and low W i = 0.01 the droplet remains nearly spherical even for the case when it is partly yielded, cementing the qualitative difference between near-viscoplastic and EVP droplets in the previous section. It is more interesting to confirm that the solid regions also increase in volume with the Weissenberg number (going from left to right in the figure). Note that a similar trend has been observed for pseudoplugs in other complex geometries 14,41 . The droplets at Bi = 0.32 and Bi = 0.8 are both fully yielded at W i = 0.01, but most of their volume is unyielded at W i = 2.31. Because our value of β = 0.68 is higher than is typical for Carbopol, we have verified that this trend persists also for β = 0.05.
Moreover to explain how this happens, we now examine the stress fields at W i = 0.5 and compare them against those at W i = 2.31 at same Bingham number (Bi = 0.8). column), first normal stress difference N 1 = τ xx − τ yy (middle column) and the deviatoric extra stress tensor |τ d | (right column). The maximum value of tr(τ ) increases slightly with W i, while the maximum value of τ d is unchanged, and for N 1 decreases slightly. However, the distribution of those maximum values confines to thinner region along the interface. Since Bi/|τ d | determines the yield threshold, when τ d decreases, the unyielded zones increase. The distributions of the tr(τ ) represent degree of local elastic deformation, which also varies more inside the droplet when W i is increased. These figures show that although the maximum stress is nearly unchanged with elasticity, the elastic droplet redistributes the stress more efficiently, leaving large central regions at lower stress and lower deformation.   We also note that in the Saramito model, the elastic stress in the unyielded regime can be written as 10,43 τ e = (µ p /λ) , where is the local deformation. The elastic stress is inversely proportional to the elastic time scale λ (i.e. W i), if the local deformation is constant. Hence, the overall stress is expected to decrease with W i, if the droplet overall deformation remains roughly constant. In our case, the droplet deformation is mainly determined by capillary forces (except in near-plastic cases at high Bi), and therefore remains similar when W i increases. Hence, the elastic stresses are indeed expected to decrease overall with W i, and the unyielded regions increase overall. As a secondary effect, it is also likely that capillary stresses at the droplet tips counteract better the external shear for the elastically deformable droplet; for the rigid near-viscoplastic droplet, the yield stress has to act as the main balancing force. To support this hypothesis, the unyielded regions are shown in Fig. 9 at different capillary numbers for droplets at two representative Weissenberg numbers (W i = 0.5, 2.31) and fixed Bi = 0.8. The capillary number increases when going from left to right in the figure, and is varied in the range of 0.1 ≤ Ca ≤ 0.5. Indeed, increasing Ca and resulting weakening of the capillary force results in an earlier yielding. A higher tip curvature is required to balance the external shear when capillary forces are weak. This implies larger deformation of the droplet, which results in higher elastic stresses than at low Ca, and therefore yielding. At W i = 0.5, the material is not elastic enough to sustain solid regions with increase in deformation with Ca, while at W i = 2.31, the droplet remains mostly solid. We also note that unyielded regions increase with W i at low Bi for all Ca considered in the current work.
A complete map of the yielded and unyielded regimes is shown in Fig. 10 for W i and Bi in the range of 0 ≤ W i ≤ 2.31 and 0 ≤ Bi ≤ 3.2, respectively. Each marker in this figure corresponds to one simulation, where the final state of the droplet has been classified into one of three regimes: i) fully yielded (hollow square), ii) partly yielded (filled square), or iii) fully unyielded (circle). The partly yielded regime denotes a droplet containing yielded and unyielded regions simultaneously. For low and moderate Bingham numbers (0 < Bi < 1.6), the droplets are in the yielded or partly yielded regimes. For high Bingham numbers however (Bi > 2.4), we find that the droplet is in the fully unyielded regime for all W i studied. For each value of W i, a critical Bingham number Bi c can be defined as the lowest Bingham where the droplet enters the partly yielded regime. The figure shows that the critical Bingham number for partly yielded regime, Bi c , decreases as a function of the Weissenberg number. This confirms that the trend observed in figure 7 is universal; the unyielded region increases in volume and penetrates further within the bulk of the drop with W i. Note that in an experiment, for a given material, an increase in the shear rate results in decreasing Bi, but an increasing W i and Re. It is possible to keep the capillary number constant experimentally by changing the droplet size. Moreover, in some experiments inertial effects are negligible. To mimic this in Fig. 10, we also included constant material curves Bi × W i = C at C = 0.1, 0.5, 1, 3 and 5, where different values of C refers to various materials. Thus following each curve from left to right implies increasing shear rate for same material, while keeping capillary number constant. On the other hand, horizontal movement between two curves implies changing the material elasticity while keeping the same shear. All droplets yield when the shear increases sufficiently, but more elastic EVP droplets yield at higher levels of shear than less elastic or near-viscoplastic droplets.

V. COMPARISON BETWEEN EVP AND PURELY VISCOELASTIC DROPLETS
While the previous section focused on the plastic and elastic effects and solid distributions of EVP droplets (information hard to find experimentally), the aim of this section is to shed light on the differences between EVP and viscoelastic droplets. The two quantities analysed here are relatively easy to measure and compare in possible future experiments: droplet deformation and orientation. This qualitative development is sharply contrasted by the VE droplet in Fig. 11(c). The deformation of the VE droplet in time exhibits a slight overshoot before saturating to a steady state value. The overshoot appears because a finite time is required to develop inhibitive viscoelastic stresses, and is therefore more pronounced for higher W i (when the characteristic time scale is longer). The steady-state value of the drop deformation changes non-monotonically with W i, which is again in agreement with the numerical results of Aggarwal and Sarkar 19 . When considering the saturated mean value of the drop deformation, in agreement with previous studies we observe that increasing W i decreases the deformation of the VE droplet, in contrast to the EVP droplet for which the deformation increases with viscoelasticity. Figure 11(b) shows the prediction of the ODE model for an EVP droplet presented in this paper (Appendix A), with a surface tension spring and an elastic spring in parallel with a viscous damper. We would like to emphasize that the model describes a fully unyielded droplet, which is the case also in the numerical simulation in (a). The model captures the development of the mean deformation with W i very well, indicating that the droplet at Bi = 8 indeed behaves like a viscoelastic solid. However, the ODE model cannot predict the time-periodic oscillations. We calculated the natural frequency of the corresponding spring-damper system, and found it to be almost an order of magnitude higher than the frequency observed in our simulations that was close to the flow shear time scale (t = 1). However, at Re = 0.05, the shear flow itself is very stable. Hence, we propose that the oscillations are a result of a fluid-structure interaction, rather than determined by either phase alone. It is worth mentioning that similar periodic oscillations were observed for viscoelastic droplets which were much more viscous than the outer fluid 21 . Here, although the viscosity ratio is low (k µ = 1.5), the droplet behaves like a very viscous one because it is fully solid.
Finally, the development of the saturated (mean) values of the drop deformation with W i is quantified in Fig. 11(d) for a wide range of Bingham numbers (Bi = 0 − 80). This figure shows how the dependence on W i changes continuously from non-monotoneous but mainly decreasing at Bi = 0, to monotoneous and highly increasing at Bi ≥ 1.6. One can observe that the viscoelastic effects are more pronounced in the unyielded regime. The long-time drop deformation at high Bi increases monotonically with the Weissenberg number and eventually reaches a plateau. The ODE prediction (orange line) is also included. Interestingly, the numerical simulation results at high Bi are close to the predictions of the ODE model, as might be expected as the model was derived for the fully unyielded regime. We cannot expect them to fully converge to each other, because of the extreme simplicity of the model.
In Fig. 12, the time development of the orientation angle θ is shown for various values of W i in the range 0 ≤ W i ≤ 2.31. The time evolution of the orientation angle is similar to the drop deformation, i.e., the EVP droplet approaches the steady state orientation significantly slower than the in viscoelastic droplet, due to overshoots accompanied with oscillations in the EVP case. Next, we will attempt to interpret the steady state values. It is known that the deformation of a Newtonian droplet increases with Ca, and that the droplet therefore aligns itself closer to the flow direction (smaller θ). The viscoelastic droplet ( Fig. 12(a)) decreases its deformation with W i, and accordingly, θ increases (further away from the flow direction). This observation is in agreement with Aggarwal and Sarkar 19 . For the EVP droplet on the other hand ( Fig. 12(b)), the trends are reversed: the drop alignment with the flow increases with W i. The change in the EVP droplet orientation leads to changes in the external flow field, and increases the flow-induced viscous stretching of the EVP droplet, as further discussed in next section.

Physical explanation of the effect of Weissenberg by normal stress differences
We now consider the distribution of stresses along the droplet interface, following the arguments of Yue et al. 18 and Aggarwal and Sarkar 19 who studied the same for viscoelastic droplets. The stress fields (Fig. 13(c)-(f)) and the drop interface ( Fig. 13(a),(b)) are shown in the central plane (z = L z /2) as functions of the angle φ measured from the streamwise direction. By analysing these distributions, we seek to understand why the trends for VE and EVP droplets are different.
The distribution of the EVP normal stress (i.e. τ pn = n · τ τ τ · n, where n is the outgoing normal of the interface) is shown for Bi = 0 (Fig. 13 c) and Bi = 1.6 ( Fig. 13 d). At both Bingham numbers they share some features. The maximum of the EVP normal stress lies in both cases in the vicinity of the drop tips while its minimum lies around the equator. Furthermore, τ pn is tensile near the drop poles, thus pulling the interface inward and reducing the long L-axis, whereas near the equator region it is compressive and pushes the interface outward, increasing the short B-axis ( figure 13c,d). Hence, the stresses at the tips and at the equator act mainly to reduce the droplet deformation.
Considering now the purely viscoelastic droplet (Bi = 0), the tensile stress at the tip shows a non-monotonic trend. In increasing the Weissenberg number from W i = 0.15 to W i = 1, the maximum value of the tensile stress increases. However, on further increasing it to W i = 2.31, the peak value slightly decreases. Because of this action, an increase in W i results in a non-monotonic trend of the drop deformation, i.e. the deformation first decreases and then slightly increases.
For the elastoviscoplastic droplet (Bi = 1.6), the tensile stress undergoes overall a much larger variation with increasing W i. The variation at the tips is still non-monotonic, but the compressive stress around the equator affecting the B-axis decreases monotoneously and significantly with W i. Because of this, the deformation of the EVP droplet increases monotoneously with W i.
Finally, the viscous normal stress is shown in Fig. 13(e,f) (i.e. τ vn = n·µ s (∇u+∇u T )·n). In general, a higher magnitude of normal viscous stress at the drop top increases L, while at the equatorial region it decreases B, both resulting in larger deformation. For the EVP droplet, it is found that the magnitude of the τ vn increases both at the drop-tip and the Summarizing what the stress distributions tell us, we observe a non-monotonic deformation due to stress distribution at the drop tips and a simultaneous increase in deformation as a result of decreasing stress at the equator. The relative magnitudes of these two phenomena cause the non-monotonic trend for Bi = 0 case and a monotonic increase for Bi = 1.6.

B. Effect of capillary number and solvent viscosity ratio
Next, we examine how the solvent-to-polymeric viscosity ratio β affects the dynamics of the droplet. In Fig. 14, the long-time saturated mean values of the drop deformation are depicted for various values of β in the range 0.1 ≤ β ≤ 0.8, while keeping the droplet-toambient viscosity ratio k µ constant. For the viscoelastic droplet in Fig. 14(a), when β is increased, the long-time value of D increases. This can be explained by that β modifies the effective Weissenberg number W i = W i(1 − β), following Aggarwal and Sarkar 20 , Izbassarov and Muradoglu 44 ; increasing β has a similar physical effect for a viscoelastic droplet as decreasing W i. Elasticity inhibits deformation for the viscoelastic droplet, hence decreasing W i increases its deformation.
Also for the EVP droplet with Bi = 1.6 ( Fig. 14(a)), on increasing β, the long-time deformation increases. By the earlier argument this seems counterintuitive, because W i decreases here thus would lead to decrease in deformation for the EVP droplet. However, it should be observed that W i is based on the first normal stress difference N 1 for the viscoelastic droplet 20 , and in Section V A 1 we could observe that the EVP droplet normal stress behaves differently. We can instead employ a simple physical argument based on the mechanistic system that describes the droplet in its unyielded state (Appendix A). The nondimensional material elasticity parameter k e = (1 − β)/W i decreases with an increase in β. Therefore for the EVP droplet, increasing β has a similar effect as increasing W i, resulting in a increase of long-time deformation, differently from the viscoelastic droplet. The effective Weissenberg number for the fully unyielded EVP droplet then becomes W i = W i/(1 − β). The proposed ODE model for the EVP droplet is also able to capture the increase in deformation (Fig. 14a). Summarizing these results, an increase in β increases the deformation for both viscoelastic and elastoviscoplastic droplets.
Finally, the effect of surface tension was investigated by varying Ca in the range of 0.1 ≤ Ca ≤ 0.5. The long-time values (Fig. 14(b)) of D are shown. In Fig. 14(b), the long-time deformation is compared with experimental data for a Boger fluid droplet in a Newtonian matrix 27 , and numerical results for an elastic particle in a Newtonian fluid system 31 . The viscoelastic droplet agrees very closely with the other two, indicating that the surface tension has the same effect on all similar viscoelastic systems (an elastic particle could be considered as a viscoelastic material with an infinite relaxation time, hence very high W i). It is worth noting that the long-time deformation of the EVP droplet becomes independent of Weissenberg number when W i is sufficiently high (Fig. 11), which may explain the good match with the elastic particle. The deformation of the EVP droplet increases similarly with Ca, but the slope is slightly smaller than for the viscoelastic droplet, and again, the effect is well reproduced by the proposed ODE model.

VI. CONCLUSIONS
A non-Newtonian elastoviscoplastic (EVP) droplet surrounded by a Newtonian simple shear flow was studied by three-dimensional direct numerical simulations. An elastoviscoplastic droplet exhibits simultaneously viscoelastic and plastic behaviours, and we studied its time-evolution under constant shear, and its steady state. First, we performed a numerical experiment with a droplet of Carbopol under increasing shear rates. The physics of the yielding process was then analysed by varying the nondimensional parameters of the problem. The yielding process in time for a near-viscoplastic EVP droplet (W i = 0.01) was much faster and qualitatively different from more elastic EVP droplets. Also, the unyielded regions (solid plugs and pseudoplugs) within the droplet increased in size with increasing Weissenberg number. This can be explained by that the elastic droplet can deform and redistribute its stresses more efficiently in the unyielded state, leaving large regions under low stress, and allowing capillary forces to act at the drop tips. As unyielded regions are larger at higher W i, this leads to a somewhat surprising conclusion: plastic effects are enhanced by droplet elasticity. A two-dimensional regime map was constructed as a function of Weissenberg and Bingham numbers, classifying each case into one of three regimes, depending on whether the droplet is i) fully yielded, ii) partly unyielded of iii) fully unyielded. In summary, our simulations predicted that the volume of the unyielded region inside the droplet increases with the Bingham number and the Weissenberg number, while it decreases with the capillary number at low Weissenberg and Bingham numbers.
Furthermore, we compared droplet deformation between viscoelastic and EVP droplets, a quantity which would be easy to measure in experiments. We found that the EVP droplet behaves in many ways differently from a viscoelastic droplet, but with a continuous transition between the two behaviours as a function of the Bingham number. Firstly, the deformation and stresses of the EVP droplet oscillated periodically in time. Such oscillations are not seen for viscoelastic droplets, but have been observed in systems with a high viscosity contrast. Secondly, the EVP droplet deformed more with an increasing Weissenberg number, in contrast with the viscoelastic droplet. The drop inclination angle also had opposite trends for VE and EVP droplets (not shown for brevity). However, both droplets deformed more with an increasing capillary number, and with a decreased polymeric viscosity (or more precisely, an increasing solvent-to-total viscosity ratio). While the dynamics of the viscoelastic droplet is governed by its first normal stress difference, the dynamics of the EVP droplet in its fully unyielded regime is mainly governed by its elastic spring constant. Finally, we proposed a simple ordinary differential equation, which captures well the observed trends of the EVP droplet parameter dependence in the fully unyielded regime. However, the model was not able to capture the time-periodic oscillations observed in our numerical simulations. In summary, while any droplet of a yield-stress fluid will yield at high enough shear (low enough Bingham number), we found that also other parameters have a significant qualitative influence on its yielding behaviour and deformation. In particular, droplet elasticity and plasticity both need to be accounted for.
In the future, we would welcome experimental studies of EVP droplets in shear flows, to confirm to which extent the qualitative differences between VE and EVP droplets can be observed in experiments. It would also be interesting to perform studies of EVP droplets with further improved models, including other physical effects that complex fluids often represent, such as shear-thinning 43 , kinematic hardening 45 or even thixotropy 46 .

VII. ACKNOWLEDGMENTS
We acknowledge financial support by the Swedish Research Council through grants No. VR2013-5789 and No. VR2017-4809. The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 852529).
Aggarwal and Sarkar 19 developed a first-order ordinary differential equation to describe the deformation of a viscoelastic drop in a steady shear flow, while here, a similar model is presented to explain the EVP drop deformation in the same flow. Eq. A1 can be rewritten in terms of forces by multiplying it with the area R * 2 as: where the spring constant k * Γ is defined as k * Γ = Γ * /R * . Eq. A2 represents a force balance between the imposed shear force from the outer flow (right-hand side), and the viscous, interfacial and EVP forces of the droplet (left-hand side). On the right hand side, F total = µ * 1 R * 2γ * term is the viscous stretching of the droplet by the outer shear flow. On the left hand side, F viscous = µ * s,2 R * 2ε * is the viscous damping force inside the droplet, F interf acial = k * Γ R * 2 ε * is the force due to the surface tension, and the last term on the left F EV P = k * e R * 2 ε * is the contribution of elastoviscoplastic stresses. In order to solve the ODE model, an initial condition is required, and an initially spherical drop ε * (0) = 0 is assumed here (as in our numerical simulations). Eq. A2 can be non-dimensionalized with the length scale R * and with the time scaleγ * −1 , thus obtaining dε dt + 1 β where Ca = µ * 1 R * γ * /Γ * , W i = λ * γ * , β = µ * s,2 /µ * 2 , k µ = µ * 2 /µ * 1 and ε are the model capillary number, the model Weissenberg number, the ratio of the solvent viscosity to the total viscosity, the ratio of droplet to bulk viscosity, and the drop deformation, respectively.
An analytical solution can be obtained from (A3), i.e.