1 Introduction

Shear dynamics of an inverted nematic emulsion

A. Tiribocchi, M. Da Re, D. Marenduzzo, E. Orlandini

Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX
First published on the web Xth XXXXXXXXXX 200X

DOI: 10.1039/b000000x

Here we study theoretically the dynamics of a 2D and a 3D isotropic droplet in a nematic liquid crystal under a shear flow. We find a large repertoire of possible nonequilibrium steady states as a function of the shear rate and of the anchoring of the nematic director field at the droplet surface. We first discuss homeotropic anchoring. For weak anchoring, we recover the typical behaviour of a sheared isotropic droplet in a binary fluid, which rotates, stretches and can be broken by the applied flow. For intermediate anchoring, new possibilities arise due to elastic effects in the nematic fluid. We find that in this regime the 2D droplet can tilt and move in the flow, or tumble incessantly at the centre of the channel. For sufficiently strong anchoring, finally, one or both of the topological defects which form close to the surface of the isotropic droplet in equilibrium detach from it and get dragged deep into the nematic state by the flow. In 3D, instead, the Saturn ring associated with normal anchoring disclination line can be deformed and shifted downstream by the flow, but remains always localized in proximity of the droplet, at least for the parameter range we explored. Tangential anchoring in 2D leads to a different dynamic response, as the boojum defects characteristic of this situation can unbind from the droplet under a weaker shear with respect to the normal anchoring case. Our results should stimulate further experiments with inverted liquid crystal emulsions under shear, as most of the predictions can be testable in principle by monitoring the evolution of liquid crystalline orientation patterns or by tracking the position and shape of the droplet over time.

footnotetext: Dipartimento di Fisica e Astronomia and Sezione INFN di Padova, Universitá di Padova, Via Marzolo 8, 35131 Padova, Italy.
SUPA and The School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, United Kingdom

## 1 Introduction

Dispersions of particles in a host fluid are important examples of soft matter with many potential applications in food industry, drugs, paints and design of new composite materials 1. While in colloidal suspensions the hosted particles are solid, emulsions are dispersions of liquid droplets coated with a surfactant. When the particles are dispersed in a nematic liquid crystal, i.e., an anisotropic fluid, where elongated organic molecules are, on average, aligned along a common direction (director), additional long-range forces, due to elastic deformations of the director field in proximity of the droplet surface, are induced and topological defects are observed. Defects mediate anisotropic droplet-droplet interactions of dipolar or quadrupolar type that typically lead to the formation of ordered droplet structures, such as chains and hexagonal lattices 2, 3, 4, 5, 6.

Due to a large body of theoretical and experimental works, the equilibrium properties of isotropic droplets in a nematic host (inverted nematic emulsions) are now well understood: depending on the strength and direction of the anchoring of the nematogens at the droplet surface, many equilibrium structures occur, each with a well defined defect and droplet shape conformation 6. For example, isotropic droplets with normal anchoring of the director field at their surface are often accompanied by a hyperbolic hedgehog defect 2, 3 while, for weaker anchoring a Saturn ring defect (or a pair of antipodal defects of topological charge -1/2 in 2D) surrounding the droplet is observed 7, 8, 9. The processes of defects formation and director field orientation in proximity of the droplet surface have been also investigated numerically by different means such as molecular dynamics simulations 10, 11, 12, 13, Monte Carlo simulations 14, 15 and through minimisation of a Landau-de Gennes free energy functional 17, 16, although in many cases the approximation of undeformable droplet has been considered.

Much less is known, however, on the rheological and response properties of isotropic-nematic mixtures when they are, for example, subject to external perturbations such as electric, magnetic and flow fields. In fact, although molecular dynamics approaches are in principle able to describe the time evolution of these systems, they are limited to very small length and short time scales where hydrodynamic modes are still irrelevant. Continuum models, based on a Landau-de Gennes free energy description of the nematogens, have been used in the past to study either the phase separation dynamics of symmetric nematic-isotropic mixtures 18 or the effect of hydrodynamic flow on the orientational order of the nematic liquid crystal in presence of spherical, rigid, inclusions 19. More recently, other theoretical studies, based on a free-energy description of the system and Lattice Boltzmann approaches, have been carried out with the aim of characterising either the equilibrium properties of the interface between nematic and isotropic fluids or the shape of a droplet of an isotropic fluid immersed in a nematic liquid crystal in presence of a surfactant 21, 20. In particular Sulaiman et al. 22 have recently introduced and tested a Lattice Boltzmann algorithm which solves numerically the hydrodynamic equations of motion of a nematic coexisting with an isotropic phase, either in absence or in presence of an electric field.

Here we adapt this algorithm to study the effect that an externally imposed shear flow can have on the dynamical properties of inverted nematic emulsions, described as a single two- or three-dimensional droplet of isotropic fluid surrounded by a nematic liquid crystal. By varying the shear rate and the ratio between the elastic energy scales of the nematic in the bulk and at the droplet surface, we observe a rich dynamical response. For weak anchoring, the behaviour of the system resembles that of a sheared isotropic droplet in a binary fluid, whereas for intermediate and strong anchoring many other steady states are observed. Of particular interest is the occurence of oscillatory steady states in which the droplet tumbles and deforms in the flow (for intermediate anchoring), or in which the topological defects created by the anchoring detach from the droplet surface and move around in the bulk. Such oscillatory steady states are possible since the system is driven far from equilibrium by the applied shear. In our work, we characterise the steady states in terms both of the droplet shape (aspect ratio and tilt) and of the defect textures within the nematic host.

The paper is structured as follows. In Section 2, we introduce the model which describes the equilibrium phase behaviour and the hydrodynamics of 2D droplets of isotropic fluid suspended in a lyotropic nematic liquid crystal. In particular we write down the Landau-de Gennes free energy of the system and the Beris-Edwards equations of motion for nematic liquid crystals coupled to Cahn-Hilliard dynamics for the diffusion of the two species. In Section 3 we partition our results as follows. First we discuss the equilibrium properties of a single droplet of Newtonian fluid dispersed into a nematic host fluid, for different anchoring conditions (either homeotropic, i.e., normal, or homogeneous, i.e., tangential). We next present the main findings of the paper referring to the effect that linear shear flow has on the shape and dynamical properties of the emulsion. The end of the section is devoted to extend the study to 3D systems where the defect configuration at equilibrium is given by a Saturn ring hugging the droplet. The effect of the shear on the dynamics of the droplet-Saturn ring pair is presented and compared with the 2D counterparts. Finally Section 4 is devoted to a discussion of the results and to conclusions.

## 2 Model and methods

### 2.1 Underlying free energy, and equations of motion

We consider an inverted nematic emulsion in an extremely diluted regime in which a single isotropic droplet is dispersed in a nematic liquid crystal. In this two-component system the droplet is made up of anisotropic-shaped molecules (such as rods) randomly distributed and oriented in the space. The resulting phase is that of a liquid crystal in the isotropic phase in which, unlike in the nematic phase, there is neither orientational nor positional order. The physics of the system can be described in terms of a set of coarse-grained variables , , and which are respectively: an order parameter related to the relative concentration of the nematic phase ( in the isotropic phase); the mass density; the average velocity field and the tensor order parameter that, within the Beris-Edwards theory 23, describes the nematic phase. More specifically, in the uniaxial approximation , where is the director field (Greek subscripts denote Cartesian coordinates) and is the local degree of nematic order related to the largest eigenvalue of ().

The equilibrium properties of the system are encoded in a Landau-de Gennes free energy , where

 f=fbf(ϕ)+flc(ϕ,Q)+fint(ϕ,Q)+fW(Q). (1)

The term

 fbf(ϕ)=a4ϕ2(ϕ−ϕ0)2+κ2|∇ϕ|2 (2)

stems from a typical binary fluid formalism and is made by two contributions: the first one is a double well potential allowing bulk phase separation into a nematic (outside with ) and isotropic (inside with ) phase in the droplet geometry, whereas the second one creates an interfacial tension between these phases whose strength depends on . The term

 flc(ϕ,Q) = A0[12(1−ζ(ϕ)3)Q2αβ−ζ(ϕ)3QαβQβγQγα (3) + ζ(ϕ)4(Q2αβ)2]+K2(∂αQβγ)2,

is the free energy density of the liquid crystal phase. It is made by four contributions (summation over repeated indexes is assumed). The first three, multiplied by the positive constant , are the bulk free energy density for an uniaxial liquid crystal system with an isotropic-nematic transiton at . The parameter , which determines which phase (isotropic or nematic) is the stable one, is assumed to be linearly dependent on the concentration , namely

 ζ=ζ0+ζsϕ, (4)

where and are constants controlling the boundary of the coexistence region. The fourth term creates an elastic penalty for local distortion of the nematic order, within the (standard) “one elastic constant” approximation 24, with being the resulting single elastic constant. The term

 fint(ϕ,Q)=L(∂αϕ)Qαβ(∂βϕ) (5)

takes into accout the anchoring of the nematic liquid crystal on the surface of the droplet. The constant controls the anchoring strength: if negative the director is aligned perpendicularly (homeotropic anchoring) to the surface, whereas if positive the director is aligned tangentially to the surface (planar anchoring). Finally, in presence of confining walls, one has to take into account the anchoring of the nematic with these boundaries. This is described by the last term of Eq. (1)

 fW(Q)=12W(Qαβ−Q0αβ)2, (6)

where the constant controls the strength of the nematic anchoring at the walls and with and being respectively the direction and the magnitude of the nematic ordering at the walls.

The dynamical equations governing the evolution of the system are

 ∂tϕ+∂α(ϕuα) = ∇(M∇δFδϕ), (7) (∂t+→u⋅∇)Q−S(W,Q) = ΓH, (8) ∇⋅→u = 0, (9) ρ(∂t+uβ∂β)uα = ∂βσtotalαβ. (10)

The first equation, that governs the time evolution of the concentration , is a convection-diffusion equation for a model B 25 where is a thermodynamic mobility parameter and is the chemical potential.

The dynamics of the liquid crystal order parameter, the tensor , is described by Eq. (8) which is a convection relaxation equation. The first two terms on the left hand side are the material derivative. Moreover, since for rod-like molecules the order parameter distribution can be rotated and stretched by flow gradients, a further contribution is needed 23. Its explicit expression is

 S(W,Q) = (ξD+ω)(Q+I/3)+(Q+I/3)(ξD−ω) (11) − 2ξ(Q+I/3)Tr(QW),

where and are the symmetric and antisymmetric parts, respectively, of the velocity gradient tensor and is the unit matrix. The constant takes into account the aspect ratio of the molecules of a given liquid crystal and determines the dynamical behaviour of the director field under shear, in particular whether it is flow aligning (), giving a stable response, or flow tumbling (), generating an unsteady response. On the right hand side is the collective rotational diffusion constant and

 H=−δFδQ+I3TrδFδQ. (12)

is the molecular field (this is the analogue for of the chemical potential ).

Finally, the last two equations are respectively the continuity and the Navier-Stokes equations for an incompressible fluid. Here is the total hydrodynamic stress which is the sum of three contributions. The first term is the viscous stress

 σviscαβ=η(∂αuβ+∂βuα), (13)

where is an isotropic shear viscosity 23. The second one is the contribution to the stress due to the liquid crystalline order and is given by

 σlcαβ = −Pδαβ−ξHαγ(Qγβ+13δγβ)−ξ(Qαγ+13δαγ)Hγβ + 2ξ(Qαβ−13δαβ)QγμHγμ+QανHνβ−HανQνβ,

where the pressure is given by (where is the temperature), and, except in proximity of the droplet surface and of the defect cores, is costant in our simulations to a very good approximation. The last term is the sum of the interfacial stress between the isotropic and the liquid crystal phases with the elastic stress due to the distortions within the liquid crystal phase

 σsαβ=−(δFδϕϕ−F)δαβ−δFδ(∂βϕ)∂αϕ−δFδ(∂βQγμ)∂αQγμ. (15)

### 2.2 Numerical aspects and mapping to physical units

Eqs. (7,8,9,10) are solved numerically by using a hybrid lattice Boltzmann method previously used in similar systems such as binary fluids 26, liquid crystals 27, 28 and active matter 29, 30, 31. Unless explicitely stated otherwise, most of the simulations are performed on a two-dimensional rectangular lattice (, ) in which an isotropic droplet is initially placed at its centre and surrounded by a nematic liquid crystal. The entire system is sandwiched between two parallel flat walls. We choose a rectangular box with a longer size along the shear direction in order to minimize the possible effects of interference between periodic images of the droplet moving with the flow.

The concentration field is initially set to zero inside the droplet and to a constant value in the bulk nematic phase. Similarly, the order parameter is initially set to zero inside the droplet and different from zero elsewhere. In particular the initial direction of the director in the nematic phase equals the one imposed at the walls. This means that the director in the bulk is along the -direction for homeotropic (or perpendicular) anchoring at walls and along the -direction if homogeneous (or tangential) anchoring is instead considered. We incidentally note that these initial conditions lead to a droplet whose final size, after equilibration, depends on the anchoring conditions of the director on its surface. A similar effect has been observed in Ref. 32 in which an isotropic droplet is embedded in a polar liquid crystal, and has been ascribed to the deviation of the equilibrium values of of the two phases from the ones at the minima of the free energy (which are ). This deviation however is, in terms of the radius of the droplet, and does not affect in an appreciable way the dynamics of the droplet under shear. Finally at the walls we impose no-slip boundary conditions for the velocity field, neutral-wetting (meaning that there are no flows of matter across the walls 26) for the concentration field and strong anchoring of the nematic (this is achieved by setting the anchoring strength ).

Starting from this initial set up, the system is let to relax into its equilibrium state that afterwards we shear by moving both walls along opposite directions at constant speed. In most simulations the following parameter values have been used: , , , , and . Note that the parameter , appearing in Eq. (8), tunes the dynamical response of the director field under shear. We have set which is a value within the flow aligning regime where, in a linear velocity profile, the director field reaches a steady state characterised by a given flow aligning angle (the Leslie’s angle) 24. For completness we have also perfomed simulations with . This corresponds to the flow tumbling regime where a steady solution under a constant shear rate no longer exists and the director rotates in a direction consistent with the vorticity of the flow. Although the overall dynamics of the director field in both regimes is diverse, we have not found appreciable differences in the defect dynamics under a shear flow and here we only report the flow aligning case.

The key parameters for determining the shape of the droplet are the tension coefficient, which we have kept fixed to (for stability reasons), the interface cross-gradient coefficient , which ranges from to (negative values for strong homeotropic anchoring and positive values for strong tangential anchoring), and the elastic constant whose values varies between and . A suitable dimensionless quantity which characterises the shape of the droplet is the ratio , where (measured in units of ), is the perimeter of the droplet (see Appendix 1 for the calculation of when the droplet shape is deformed under shear) and its radius. This quantity measures the strength of the surface anchoring relative to the bulk elastic deformations 24. Typical values for the anchoring strength found in the literature 33, 34 range from J m to J m, for an elastic constant of N and for a droplet size of -m. This leads to a corresponding value of varying within the interval -. If surface anchoring is weak compared to bulk nematic distortions, whereas if it is larger and can significavely deform the director near the droplet. A further quantity which affects the droplet shape is the surface tension , where . If , in simulation units. This corresponds to a physical value of J m, whereas experimental values vary in the range J m 35, 33, 34. Unless stated otherwise, we have kept this quantity constant. In the next section we show that, depending on the values of these parameters, several possible equilibrium shapes, whose structure is affected by the position of topological defects, can be identified. These equilibrium states will in turn respond differently when sheared. Similar shapes have been found and discussed in Ref. 22.

Note that all the aforementioned values are in simulation units. In order to map them into physical units we follow the approach given in Ref. 22. In particular, assuming to model a flow-aligning regime, m thick device with a rotational viscosity of roughly poise (a typical value of 5CB), one gets the length-scale and time-scale to be respectively m and s. Furthermore, an elastic constant ranging between in simulation units, corresponds to pN, within the typical values of nematic liquid crystals. To compare the effect of anchoring in simulations and experiments, we can use the dimensionless parameter defined previously.

## 3 Results and discussion

We first characterize the equilibrium properties of an isotropic droplet in nematic when either homogeneous or homeotropic anchoring of the director at the surface of the droplet is considered. For each case we also look at the effect that different directions of the anchoring at the walls (i.e. either perpendicular or parallel) may have on the equilibrium configurations.

Later on the dynamical response of the equilibrated configurations when subject to a shear flow (flow aligning regime) is studied for different shear rates.

### 3.1 Isotropic droplet in a nematic phase at equilibrium

We first consider an isotropic droplet of radius located at the centre of a rectangular box of size , . The anchoring on the surface of the droplet is homeotropic and the director is perpendicularly anchored on both walls. The corresponding simulation is run for time-steps until the droplet and the surrounding medium are completely equilibrated, a state achieved when the total free energy is at its minimum. In Fig. 1 we show four equilibrium configurations obtained by appropriately changing the elastic constant and the surface anchoring strength . In terms of the number , these correspond to (Fig. 1a), (Fig. 1b), (Fig. 1c), (Fig. 1d)***Since is negative we will implicitly consider its absolute value..

When strong homeotropic anchoring is set on the surface of the droplet (see Fig. 1), an imaginary defect of topological charge is enucleated at its centre. The conservation of the topological charge requires though the formation of two defects of topological charge , which are located on opposite parts of the droplet (see, for instance, Fig. 1c or d). This is the 2D version of the well-known Saturn ring 2, 3, 36Other configurations are possible, such as a hyperbolic hedgehog (a defect of charge located on one side of the droplet) or a disclination ring of a finite radius (located above or below the droplet). An extensive discussion can be found in Ref. 36. The position of the defects pair can be controlled by properly balancing the strength of the surface anchoring and the bulk elastic distortions (namely the parameter ). For , for instance, surface anchoring is very weak and defects disappear, leaving the shape of the droplet unaltered (Fig. 1a). Notice that this is in agreement with the topological charge conservation, as in absence of anchoring the imaginary defect (of charge ) does not form inside the droplet, maintaining the total charge zero. When the anchoring strength becomes comparable with the elastic nematic energy (), defects appear on the surface where they are firmly anchored (Fig. 1b), favouring a nutshell-like shape of the droplet. For higher values of , both defects emerge on opposite sides of the droplet along the equator (Fig. 1c). Lastly, for , since the elastic liquid crystal energy becomes very small compared to the surface energy, the defects pair are clearly far apart from the droplet and well inside the bulk nematic phase. Similar droplet configurations have been also discussed in Ref. 22 in which their dynamics in presence of an applied electric field is studied, and in Ref. 20 where the role of a surfactant has also been taken into account.

For homogeneous anchoring at the droplet surface, different steady states are expected. In this case we have considered an isotropic droplet of radius embedded in the nematic phase with homogeneous anchoring on both walls. This choice determines a final configuration in which two defects of topological charge are located on opposite sides along the equator, similarly to what obtained for the states of Fig. 1. By varying the values of the nematic elastic constant we have identified two equilibrium shapes (see Fig. 2). When the distortion energy overcomes the surface one, defects are completely absorbed within the droplet and the director profile smoothly surrounds its surface (Fig. 2a). When both effects become comparable, two defects form on both side of the droplet, keeping unaltered the total topological charge (Fig. 2b). The position of defects changes if homeotropic anchoring is set on both walls. In this case they will form along the north-south direction of the system, on opposite sides of the droplet (not shown). However this final state is unstable to small perturbations and both defects eventually shift and move along the droplet surface until they find a more stable conformation.

We finally mention that, in line with previous nematodynamics studies 37, 38, we found a small degree of biaxiality especially close to defects. By following the approach used in Ref. 38, regions of biaxial order can be found by calculating the values of three parameters, namely , and , where , and (with ) are the eigenvalues of the diagonalised matrix . These parameters have the following properties: and . A well-ordered uniaxial nematic arrangement will give , whereas regions of isotropy and of planar ordering (the biaxial state) correspond to and , respectively. In particular we found near the defects and far from them, either in the nematic or in the isotropic phase.

### 3.2 Isotropic droplet in a nematic phase under shear

By starting from the equilibrated droplets previously described, we now impose a shear flow on the system by moving the top wall along the -axis with velocity and the bottom wall in the opposite direction with velocity . This sets a shear rate measured in in simulation units. In all the cases studied, for low shear rates (typically for ) the initial droplet configuration is the same as for an unsheared system, independently on the values of the elastic constant and on the surface anchoring constant . We will therefore discuss only results whose shear rate is strong enough to induce deformations and/or droplet motion.

Besides the parameter , the shape of the droplet and the role played by the forces in the system can be characterized by introducing the capillary number . This adimensional quantity, often used in rheological experiments, measures the strength of viscous forces relative to the surface tension acting at the interface between two immiscible fluids. In these simulations is expected to range from to . To quantify the effect of the shear on the droplet shape we also consider the parameter  39 ( and represent the major and the minor axis respectively), which measures the droplet deformation under shear and goes from (no deformation) to (full deformation). A further parameter worth considering is the Reynolds number , which measures the importance of inertial forces relative to the viscous ones: for low shear but, as the shear rate increases, increases up to where inertial forces become comparable with the viscous ones and the condition of laminar flow is less complied with.

A well-known result in binary fluid dynamics is that, if an immiscible isotropic emulsion is dispersed in a Newtonian fluid, one would expect that, under shear, the droplet deforms and orients along the shear flow. Therefore, a first check to validate our model consists in simulating a sheared system in which an isotropic droplet with no surface anchoring (, or ) is embedded into a nematic liquid crystal having a modest or low elastic energy. This is the closest approximation of our system to the simple binary fluid case. In Fig. 3a we show the steady state attained by the droplet after imposing a shear (with ) and the director profile of the nematic phase. Similarly to what observed in immiscible mixtures of Newtonian fluids, the droplet elongates and afterwards aligns along the direction of the shear flow, with a deformation rate that, at steady state, is . The velocity field (shown in Fig. 3b) is more intense near the walls and weak in the central region of the lattice, with a clockwise vortex generated by the droplet. The angle formed by the major axis of the droplet with the shear direction is degrees, very close to values reported in literature for a single isotropic phase droplet under shear at low Reynolds number (for instance in Ref. 41 the authors found an angle of degrees for ). Note that in our case and . These numbers indicate that inertial forces are still weak and the whole system operates in the Stokes regime. In addition, the director field aligns along the shear flow far from the droplet and in its neighbourhood follows almost everywhere the direction of the major axis. These results will be assumed as benchmark to be compared with the other cases discussed in the following sections, in which the surface anchoring is no longer neglected (meaning ) and higher shear rates are considered.

Finally, we mention a couple of subtle points which need to be kept in mind in setting and interpreting our simulations. First, it is well known that, by increasing the shear-rate, the temperature at which the isotropic-nematic transition occurs decreases 40. This would correspond to a conversion of the isotropic phase of the droplet into a nematic one. In our simulations we avoid this effect by keeping the shear rate small enough. Second, one needs to avoid interface-interface interactions that can occur in strongly deformed droplets under strong shear. This is achieved by considering a droplet with a sufficiently large radius (provided that its interaction with walls at equilibrium remains negligible) embedded in a relatively long rectangular lattice (necessary to diminish the effect of the periodic image of the droplet). In our simulations, as the typical interface thickness separating the droplet and the liquid crystal is around lattice sites (which corresponds to roughly m), a radius of at least lattice sites is needed. In the next Sections we will show the results obtained for droplets with a larger radius () and compare these with the ones observed for a smaller radius (). In the former, the equilibrated droplet states are the same as the ones seen for the smaller radius even though obtained for slightly different values of .

#### 3.2.1 Weak homeotropic anchoring

We initially consider the case in which homeotropic anchoring is set on the surface of a droplet. This anchoring can be experimentally achieved by means of chemical treatments, such as coating the droplet with a surfactant which favours a perpendicular alignment. As previously mentioned, we have simulated a droplet of radius located in a rectangular box of size , (area fraction at equilibrium ) and a droplet of radius in a box of size , (area fraction at equilibrium ).

When the equilibrated droplet assumes the shape of Fig. 1a (independently of the radius). The dynamics under a moderate shear (namely for ) is expected to be very similar to the case discussed in Fig. 3, in which a droplet with no surface anchoring () has been studied. Indeed for , besides the different elastic constant (now ) and the weak contribution of the surface anchoring (), the shear stretches and elongates the droplet along the shear flow without generating any net motion We set necessary to avoid an excessive shrinkage of the isotropic phase (droplet), hence at the steady state.. In particular the droplet reaches a steady state with the major axis forming an angle degrees with the shear direction. The steady state values of the deformation and capillary number are respectively and , not very far from the ones obtained for the case. A very similar dynamical behaviour is observed for . In this case an angle degrees is attained at the steady state when , with a deformation rate and a capillary number

An increase of the shear rate determines a shrinkage of the droplet (regardless of its size) and may lead to its disappearance for very high . However larger droplets permit the study of the physics in a wider dynamical range and unveil unexpected properties. For instance, a well-know effect observed on a droplet in a binary fluid mixture is its break-up as the shear is increased 42. This phenomenon is still observed when isotropic droplets are surrounded by a nematic liquid crystal, regardless of the presence of topological defects. In particular, we have identified two possible rupture regimes related to two different values of . For the droplet breaks in two smaller droplets whereas for higher values, for instance , three smaller droplets result. In Fig. 4 we show the dynamics of the former case. The droplet, initially moderately stretched along the shear direction, afterwards undergoes a deeper deformation (see Fig. 4a-b) which further elongates it. Due to the high shear rate, the rod-like state exhibited in Fig. 4b is however unstable, and the droplet breaks (Fig. 4c) and divides in two separate smaller droplets (Fig. 4d) which move along opposite directions. Before the rupture the rate of deformation increases up to . Unlike the single droplet case in which the velocity field displays one vortex inside the droplet (as in Fig. 3b), now two vortices are formed in the two smaller droplets, with a magnitude much lower than the one near the walls. At late times both resulting droplets gradually shrink and completely disappear.

#### 3.2.2 Intermediate and strong homeotropic anchoring: the role of defects

Due to the low value of the anchoring strength , the dynamics discussed so far shows several similarities with the one observed in a single isotropic droplet (binary fluid-like) under shear. On the other hand, when the droplet attains a nutshell-like steady state in which two defects are located on its surface (see Fig. 1b). This significatively affects its dynamics under shear, as shown in Fig. 5 and in Movie S1. If , the droplet, similarly to the previous case, stretches and elongates along the shear direction (Fig. 5a-b). The tilt of the major axis of the droplet can be measured by looking at the time evolution of the angle it forms with the shear direction, as reported in Fig. 6b (left scale): the droplet orientation angle initially achieves a maximum at degrees and afterwards relaxes to a constant value of degrees. Notice that the value of the angle at the steady state is very similar to the one measured for the case where defects are absent. During their rotational relaxation the two defects migrate in phase along the surface, but out of phase with the rotational motion of the major axis of the droplet, before getting pinned on opposite sides at the steady state. This difference might be ascribed to the flow generated by the symmetric shear. This stretches the droplet and changes the local orientation of the director field near the defects, pushing them along opposite directions. While the axis connecting them rotates of an angle larger than degrees around the centre of mass of the droplet, that of the droplet aligns along the flow direction, almost parallel to the director field in the bulk. It is worth noting that defects are located on the surface of the droplet but not belonging to it as they are part of the liquid crystal phase. In addition the droplet acquires unidirectional motion along the -axis with an almost constant speed (see Fig. 5c-d) as the position of the component of the centre of mass is slightly shifted upwards (see Fig. 6a in which the and the components of the centre of mass are reported), and attains a final steady state whose shape is only weakly deformed by the shear (Fig. 6b, green crosses, plots the deformation parameter ). We call it bound state (BS), as both defects remain on the droplet surface.

A persistent rotation of the defects can be achieved if the shear rate is further increased. Movie S2 shows the dynamics with . Similarly to the previous case, the droplet quickly deforms and aligns along the shear direction while both defects rotate clockwise. However, due to high value of , the velocity field (very intense near the walls but weak in the centre, where the typical vortex pattern inside the droplet forms) pushes the defects over the position attained for lower and drives them around the entire surface of the droplet. In particular, during an entire cycle, they speed up their motion near the walls of the cell and then, at the end of each cycle, slow down as their respective distance from the walls augments. This dynamics also leaves temporary wake-like signatures on the director field departing from the defects (red wakes in Movie S2) and short-living opposite charge topogical defects which annihilate quickly (red spots appearing in the nematic phase, see Movie S2). A measure of the angle that the major axis forms with the shear direction is reported in Fig. 7 with a plot of the elastic free-energy

 Fel=K2(∂αQβγ)2+L(∂αϕ)Qαβ(∂βϕ). (16)

At long times the angle displays two close local maxima, stabilized around degrees, spaced out by two local minima, one short and one large, both around degrees. Although defects continuously rotate on the surface of the droplet the major axis exhibits a characteristic oscillatory behaviour, in which the two minima are achieved when defects slow down their motion (far from the walls) and the two close maxima during the successive cycle (near the walls) (see Fig. 7, red plusses, left scale). The elastic free energy of the nematic oscillates as well at long times (see Fig. 7, green crosses, right scale); in particular its local maximum corresponds roughly to the large minimum of and its minimum to the local maximum (on its left) of observed during a cycle. A crude explanation of this can be arguably related to the backflow: the velocity field, much higher near the walls than in the middle of the cell, aligns the director field more strongly when defects are closer to the walls (hence diminishing the elastic free energy) then far from them. Indeed, the director profile observed in these states supports this interpretation (see the insets (c) and (d) of Fig. 7). On the other hand, the director field at the other extremes of the angle is significatively different (see the insets (a) and (b) of Fig. 7). In particular the state (b) has a lower elastic free energy than that in (d), meaning that elastic deformations are weaker in (b) than in (d) (notice that these states correspond to the two minima of ) whereas the elastic free energy of the state (a) has a similar value of that in (b). Interestingly, both the director and the corresponding velocity field (in the state (a)) acquire an out of plane component (along the -direction) which accompanies the escape of the cores of the defect pair into the third dimension. The out of plane component of the velocity field in particular is roughly one order of magnitude lower than the other two components. Unlike a droplet-free flow aligning nematic liquid crystal, the sole presence of topological defects on the droplet surface unveils an unexpected flow-tumbling-like dynamics which is overall akin to that observed in a strictly two-dimensional system but for a slightly higher value of the shear rate. In the latter in particular both the director and the velocity field are not significatively different from those observed in the quasi-2d case, and the final steady state (the oscillatory bound state) is preserved. In addition, at the steady state the rate of deformation of the droplet weakly oscillates around and the capillary number is . The steady state described above, in which the shear stress induces an oscillatory-rotational motion of the antipodal defects pair close to the droplet, is, to our knowledge, a new result for an inverted nematic emulsion and we call it the oscillatory bound state (OBS).

When the defects of the equilibrium configuration are further apart from the droplet surface (see Fig. 1c-d): this couples to the imposed shear to give rise to further novel dynamical behaviours. In particular, depending on the value of the shear rate, one can identify three additional dynamical regimes: for very low shear rates (tipically ) both defects remain close to the droplet surface, as in the oscillatory bound state; for intermediate shear rate only one defect goes away from the droplet (single bound (SB) state) and finally, for sufficiently high shear rates, both defects move into the bulk §§§For extremely high shear rates the isotropic phase disappears as usual. (unbound (U) state). In Fig. 8 and in Movie S3 we report the dynamics of the SB state. For the droplet initially behaves similarly to what seen for the BS case where both defects, located on opposite sides along the equator, rotate simultaneously in the clockwise direction (Fig. 8a-b). Later on the rotation is arrested and the droplet-defect system, dragged by the fluid flow, slowly moves unidirectionally along the -axis (Fig. 8c), whereas the defect near the bottom wall gradually leaves the droplet and moves into the bulk along the opposite direction (from right to left in the bottom half part of the system) (Fig. 8d). This defect motion leaves a characteristic comet-like signature of the director field in the bulk. Note that the SB state of Fig. 8d occurs because, at equilibrium (i.e. before the shear is switched on), the droplet centre of mass is not equidistant from the two walls ( as shown in Fig. 9a). This is the reason why, under shear, the droplet is slowly dragged by the flow along the -direction and, more importantly, reduces the distance between its surface and the top defect, increasing at the same time the distance from the bottom defect, which is then pushed backwards by the fluid. It is interesting to note that, after the bottom defect moves into the bulk, the droplet accelerates attaining a novel steady state with a higher velocity. During this process the total and the elastic free energies show a very similar behaviour: at early times (when defects just rotate around the droplet) they increase rapidly. At intermediate times instead they increase very slowly (this regime refers to the situation in which the defects move slowly along opposite direction parallel to the walls). Finally, at longer times, when the droplet and the bottom defect move apart, they increase quite rapidly again (see Fig. 9b).

A state in which both defects migrate into the bulk (U, unbound, state) is achieved by further increasing the shear rate. In Fig. 10 and in Movie S4 we show this situation (). The early times dynamics is similar to that observed for the SB state (see Fig. 10a-b), except for a larger deformation of the droplet ( here while in the SB state of Fig. 8). However, the shear rate is now intense enough along the -direction to detach both defects from the droplet. Interestingly, there is not an appreciable motion of the droplet along the shear direction until the defects, reappearing at the periodic side of the lattice, reapproach the droplet again. It is worth noticing that an unbound state could be also achieved by suitably controlling the elasticity of the liquid crystal. Indeed, we have found a very similar dynamics for a droplet with a smaller radius () with the parameters of Fig. 1d (hence with ) but with .

### 3.3 Diagram of non equilibrium steady states

The results just presented can be summarised and rationalised into a phase diagram in which the observed steady states are reported in the plane. The diagram displays different dynamical responses of the droplet under a linear shear. For , when defects are not present, a dynamical response typical of a Netwonian emulsion dominates for all values of shear rates . In particular if the droplet deforms and elongates along the shear flow whereas, for higher values of , it breaks into smaller droplets and eventually evaporatesWe borrow this term from the liquid-gas thermodynamics as the shrinkage of the isotropic phase is due to the shift of the temperature of the isotropic-nematic transition which favours the liquid-crystal phase.. If, on the other hand, , defects pair forms very close to the droplet surface and, due to the interplay between their dynamics and the deformation-rotation dynamics of droplet, three new interesting steady states can be identified. For small shear rates defects are bound in proximity of the surface of the droplet and no appreciable difference with the standard binary fluid behaviour is observed. Increasing unveils a region in which the motion of the droplet along the shear direction combines with a partial rotation of the defects along its surface. If is further increased, a second type of steady state occurs, in which defects are still bound to the droplet surface but now persistently rotate around it. Clearly, for sufficiently high values of , the droplet breaks up and eventually evaporates. For the defects pair forms still in proximity of the droplet but in the nematic phase and this gives rise to two additional steady states: in the first one, for intermediate values of , only one defect remains close to the droplet while the other one starts to move freely with the shear flow (single bound state). When, on the other hand, is sufficiently high, the second defect leaves the droplet as well and moves in the bulk (unbound state).

### 3.4 Homogeneous anchoring

When homogeneous (or tangential) anchoring is set on the surface of the droplet a different dynamical response of the system is observed. For these cases, our simulations are performed by using a rectangular box of size , with a droplet of radius . When defects are imaginary (see Fig. 2a), the dynamical behaviour under shear is similar to the one observed with homeotropic anchoring at low shear rate (): the droplet aligns with shear flow with . In the presence of defects at the droplet surface (as in Fig. 2b), however, the situation is different from that with homeotropic anchoring. In this case, we observe a slow drift of the droplet along the -direction induced by the defect at its right. Accordingly, the droplet increases its distance from the defect located on its left, which gradually leaves the surface of the droplet (see Movie S5). This establishes a transient single-bound state which lasts until the “free” defect wraps the periodic boundaries: when this happens, the mobile defect rotates clockwise around the surface and this time remains bound, trapped by the elastic interactions. A similar dynamics has also been observed for higher shear rates, before the droplet evaporates. Steady unbound states with defects in the bulk of the nematic host could instead be found either by decreasing the nematic elastic constant or by increasing the surface anchoring . Indeed, when imposing the equilibrated droplet looks very similar to the one of Fig. 2; however, under a shear rate of we now observe defect detachment. This is possible because a large value of augments the distance between each defect, which favours detachment.

### 3.5 Extension to 3D systems

The rich phenomenology found in the 2D system may be observed in thin film of inverted nematic emulsions under shear; however it is of interest to ask to which extent this also occurs in a fully 3D system. For simplicity, we only consider here the case where without shear the droplet is accompanied by a Saturn ring defect (the case of a hyperbolic hedgehog defect may also be of interest). Previous experimental and numerical studies on micron-size droplets hosted in a liquid crytstal and advected by a flow indicate that the Saturn ring defects formed around them are visibly displaced in the downstream direction and eventually collapse into a hyperbolic point defect 16, 43. Here we look at how a linear shear may impact on the Saturn ring-droplet system by restricting ourselves to 3D isotropic droplets hosted in a liquid crystal fluid that is homeotropically anchored to the droplet surface. The system is made by a single droplet of radius sandiwched between two walls at and . Along the and directions, periodic boundary conditions are considered with and . Since, as a first approximation, the droplet is expected to assume a generic ellipsoidal shape, we compute a set of three dimensionless numbers, , and , measuring respectively the degree of prolateness, sphericity and oblateness of the ellipsoid.

These quantities are defined as

 cl = a−ba+b+c, (17) cs = 3ca+b+c, (18) cp = 2(b−c)a+b+c, (19)

where , and are the principal axes of the ellipsoid, such that . For example, to establish whether the ellipsoid is more prolate than oblate one can compare with : if , the ellipsoid tends to be prolate whereas if the droplet is more oblate. Two extreme cases can be identified: for the ellipsoid degenerates into a segment while, if (i.e. and ) the ellipsoid degenerates into a two-dimensional circle. Finally if the spherical shape is recovered ().

Since, similarly to the 2D case, we expect that the equilibrium location of the defect and its dynamics under shear flow depends on the adimensional ratio , we consider the two extreme cases of and . In both cases the system is first let to equilibrate in absence of shear; the resulting configuration is then used as initial condition of the shear experiment.

The equilibrium state for is shown in Fig. 12a: it is readily seen that the droplet assumes a slightly deformed spherical shape with , and . The nutshell-like shape of the droplet is due to the location of the Saturn ring defect right at its surface, a situation similar to that reported in Fig. 1b for the corresponding 2D case. Note that the director field is almost everywhere parallel to the -axis, except in proximity of the droplet surface where, due to the strong homeotropic anchoring, a large splay-bend distortion emerges (Fig. 12a). When a moderate shear flow is imposed on the system (), the droplet does not move with the flow but simply stretches and elongates along the shear flow (resembling the 2D case of Fig. 5), and achieves a final steady state whose shape is that of a prolate ellipsoid (), with its major axis oriented at roughly degrees with the shear direction (Fig. 12b and Movie S6).

As expected, the director field orients preferentially along the shear direction at the centre of the system while it strongly bends both in proximity of the Saturn ring and towards the walls. The Saturn ring itself does not experience an appreciable deformation but rotates in the plane with its symmetry axis (passing through the centre of the ring) remaining almost parallel to the major axis of the ellipsoid.

For higher values of the shear rate the dynamics of the droplet and accompanying Saturn ring is significatively different. In Fig. 13 we report a simulation in which (see also Movie S7). After the shear is switched on, the droplet initially aligns and elongates along the flow direction (Fig. 13b) achieving a highly prolate () intermediate (non-steady) state (Fig. 13c). In this state the Saturn ring remains firmly anchored at the droplet surface. Interestingly, though, instead of surrounding the equator of the droplet (as for the case when ), the disclination ring follows the droplet deformation and stretches along the entire surface. Later on a more complex rearrangement is observed: the droplet moves slightly upwards (along the -axis), in regions of the system in which the flow (direct along the -axis) is more intense, and rotates around its major axis (Fig. 13d). Finally it is pushed forwards and the Saturn ring slips downstream, opposite to the direction of motion (Fig. 13e-f). The downstream motion of the Saturn ring agrees with the previous studies on similar systems 16, 43 and can be considered as the 3D concounterpart of the 2D bound state observed before (see Fig. 5).

For the equilibrium configuration is characterised by a Saturn ring surrounding the droplet but now fully embedded in the nematic phase (see Fig. 14a). As for we can distinguish two dynamic regimes, depending on the value of the shear rate. While for a moderate shear rate the Saturn ring mantains its shape unaltered and simply orients its axis along the major axis of the deformed droplet, for higher values a more complex behaviour is observed (see Fig. 14 and Movie S8, in which ). Initially the droplet strongly deforms (, and ) with its major axis aligning almost parallel to the shear flow, while the Saturn ring, although located almost at the centre of the droplet, undergoes an -like deformation (Fig. 14b-c), less pronounced though than that observed for and . Afterwards, the droplet is advected by the flow along the -axis while the disclination ring slips towards the rear part of the emulsion where eventually it gets pinned (Fig. 14d-e-f). Interestingly the shift of the Saturn ring occurs whenever the droplet acquires motion, namely when, for instance, its centre of mass shifts (either upwards or downwards). Thus, unlike in 2d, where for by increasing the shear rate, the SB and the U state could be observed, here the defect ring sticks to the droplet even for intense shear rates. The further dimension therefore appears to diminish the number of possible dynamical states the system can explore.

Note that, although the starting position of the centre of mass of the droplet is placed at the centre of the simulation box, especially for high values of the shear rate this situation is highly unstable, and any small variations either in the director field or in the velocity field inevitably drags the droplet away from it. For higher values of the shear rate we may expect the Saturn ring to further shrink into a hyperbolic hedgehog, although we have not found evidence of this effect within the simulation times considered here.

## 4 Conclusions

In conclusion, we have studied the dynamics of a 2D isotropic droplet immersed in a nematic host under an imposed symmetric shear flow. The physics of the steady states is mostly determined by the shear rate and the strength of the anchoring of the director field at the droplet surface. For weak or no anchoring, the system behaves similarly to a droplet in a binary mixture of isotropic, Newtonian, fluids. On the other hand, for intermediate or strong anchoring a variety of non-trivial non-equilibrium steady states are possible. One option is to create an oscillatory steady state, where the droplet, for instance, tumbles in the flow. Another option is to create highly dynamic states where one or both the defects, which accompany the droplet when quiescent, detach from the surface of the droplet and become mobile due to the imposed shear. For 3D systems and homeotropic anchoring the defect pair becomes a Saturn ring surrounding the droplet. We have shown that, regardless of the shear rate and of the surface anchoring strength, the disclination ring always remains localized around the droplet although it can be highly deformed and shifted downstream (similarly to the moving bound state observed in 2D). In this respect the unbound and single defect states observed in 2D can be seen as a limiting case of an inverted emulsion confined within a very thin film of liquid crystal.

All these findings suggest that the rheological response of an isotropic droplet in a nematic host is much richer than previously probed, and can further be tuned by varying strength and nature of the anchoring, together with the magnitude of the imposed flow. A viable route to test our results experimentally is, for instance, in microfluidic and rheology experiments on water-liquid crystal emulsions (or other mixtures contaning a liquid crystal and a Newtonian fluid).

## 5 Appendix 1

Here we briefly describe how we calculate the perimeter of the droplet when its shape is deformed under shear. More specifically, our simulations show that the droplet aligns along the direction of the shear flow, elongates and often acquires motion, with a shape reminding that of an ellipse. In order to calculate the ellipse which best reproduces the boundary of the droplet, we compute the inertia tensor of a flat elliptical disk, which is given by

where and are the coordinates of the centre of mass of the droplet. The mass of the -th lattice site is related to the local concentration as follows:

 mi=1−(ϕi−ϕmin)ϕmax−ϕmin, (20)

where and are the minimum and the maximum values of in the system. In particular in the nematic phase and in the isotropic phase. With this definition, in the isotropic phase and in the liquid crystal phase. In addition, for numerical stability we have imposed a cutoff on the mass values, namely we have set when . By diagonalizing the matrix , one gets the principal axes of inertia, which coincide with the and axes in an unsheared system. The principal moments of inertia are and , where and are the eigenvalues of I, , and and are the semiaxes. An estimation of the perimeter can be thus obtained through the Ramanujan formula 44

 Σ=π[(a+b)+3(a−b)210(a+b)+√a2+14ab+b2]. (21)

Lastly, the angle that the major axis forms with the direction of the shear flow can be calculated from the eigenvectors of the matrix.

## References

• 1 I. Musevic, M. Skarabot, M. Humar, J. Phys.: Condens. Matter, 23, 284112 (2011).
• 2 P. Poulin, H. Stark, T. Lubensky, and D. Weitz, Science 275, 1770 (1997).
• 3 P. Poulin, D. A. Weitz, Phys. Rev. E 57, 1 (1998).
• 4 M. Yada, J. Yamamoto, and H. Yokoyama, Phys. Rev. Lett. 92, 185501 (2004).
• 5 I. I. Smalyukh, O. D. Lavrentovich, A. N. Kuzmin, A. V. Kachynski, and P. N. Prasad, Phys. Rev. Lett. 95, 157801 (2005).
• 6 V. G. Nazarenko, A. B. Nych, and B. I. Lev, Phys. Rev. Lett. 87, 075504 (2001).
• 7 O. Mondain-Monval, J. C. Dedieu, T. Gulik-Krzywicki, and P. Poulin, Eur. Phys. J. B 12, 167 (1999).
• 8 Y. Gu and N. L. Abbott, Phys. Rev. Lett. 85, 4719 (2000).
• 9 J. Fukuda and H. Yokoyama, Eur. Phys. J. E 4, 389 (2001).
• 10 D. Andrienko, G. Germano, and M. P. Allen, Phys. Rev. E 63, 041701 (2001).
• 11 M. Allen, Comput. Phys. Commun. 169, 433 (2005).
• 12 J. L. Billeter and R. A. Pelcovits, Phys. Rev. E 62, 711 (2000).
• 13 B. T. Gettelfinger, J. A. Moreno-Razo, G. M. Koenig Jr., J. P. H.-Ortiz, N. L. Abbott, J. J. de Pablo, Soft Matter 6, 896 (2010).
• 14 S. Grollau, E. B. Kim, O. Guzman, N. Abbott, and J. J. de Pablo, J. Chem. Phys. 119, 2444 2003;
• 15 R. W. Ruhwandl and E. M. Terentjev, Phys. Rev. E 56, 5561 (1997).
• 16 C. Zhou, P. Yoe and J.J. Feng, J. Fluid. Mech. 593, 385 (2007).
• 17 H. Stark, Eur. Phys. J. E 10, 311 (1999).
• 18 T. Araki and H. Tanaka Phys. Rev. Lett. 93, 015702 (2004).
• 19 J. Fukuda, H. Stark, M. Yoneya, and H. Yokoyama, J. Phys.: Condens. Matter 16, S1957 (2004).
• 20 S. V. Lishchuk, C. M. Care, Phys. Rev. E 70, 011702 (2004).
• 21 C. M. Care, I. Halliday, K. Good, S. V. Lishchuk, Phys. Rev. E 67, 061703 (2003).
• 22 N. Sulaiman, D. Marenduzzo and J. M. Yeomans, Phys. Rev. E 74, 041798 (2006).
• 23 A. N. Beris, B, J. Edwards, Thermodynamics of Flowing Systems, Oxford University Press, Oxford (1994).
• 24 P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, 2nd Ed., Clarendon Press, Oxford (1993).
• 25 P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
• 26 A. Tiribocchi, N. Stella, G. Gonnella and A. Lamura, Phys. Rev. E 80, 026701 (2009).
• 27 A. Tiribocchi, O. Henrich, J. S. Lintuvuori and D. Marenduzzo, Soft Matter 10, 4580 (2014).
• 28 A. Tiribocchi, G. Gonnella, D. Marenduzzo, E. Orlandini, F. Salvadore, Phys. Rev. Lett., 2011, 107, 237803.
• 29 M. E. Cates, O. Henrich, D. Marenduzzo, K. Stratford, Soft Matter 5, 3791 (2009).
• 30 E. Tjhung, D. Marenduzzo, M. E. Cates, Proc. Natl. Acad. Sci. USA 109, 12381-12386 (2012).
• 31 E. Tjhung, A. Tiribocchi, D. Marenduzzo, M. E. Cates, Nature Comm. 6, 5420 (2015).
• 32 G. De Magistris, A. Tiribocchi, C. A. Whitfield, R. J. Hawkins, M. E. Cates and D. Marenduzzo, Soft Matter 10, 7826 (2014).
• 33 J. Cognard, Mol. Cryst. and Liq. Cryst. Supp. 1, London (1982).
• 34 V. J. Anderson, E. M. Terentjev, S. P. Meeker, J. Crain and W. C. K. Poon, Eur. Phys. J. E 4, 11 (2001).
• 35 S. Faetti, V. Palleschi, Journ. Chem Phys. 81, 6254 (1981).
• 36 T. C. Lubensky, D. Pettey, N. Currier, and H. Stark, Phys. Rev. E 57, 610 (1998).
• 37 N. Schopohl, T. J. Sluckin, Phys. Rev. Lett. 60, 755 (1988).
• 38 A. C. Callan Jones, R. A. Pelcovits, V. A. Slavin, S. Zhang, D. H. Laidlaw, G. B. Loriot, Phys. Rev. E 74, 061701 (2006).
• 39 G. I. Taylor, Proc. R. Soc. Lond. Ser. A, 146, 501 (1934).
• 40 P. D. Olmsted, P. M. Goldbart, Phys. Rev. A 46, 8 (1992).
• 41 Y. Y. Renardy, V. Cristini, Phys. of Fluids 13, 1 (2001).
• 42 A. J. Wagner and J. M. Yeomans Int. Journ. Mod. Phys. C 8, 773 (1997).
• 43 S. Khullar, C. Zhou and J.J. Feng, Phy. Rev. Lett. 99, 237802 (2007).
• 44 M. B. Villarino, Journ. Ineq. Pure and Appl. Math. 7, 1 (2006).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters