Morphology and flow patterns in highly asymmetric active emulsions

TL;DR

This study uses numerical simulations to explore how activity influences the morphology and dynamics of highly asymmetric active emulsions, revealing transitions from ordered droplets to complex patterns depending on activity type and level.

Contribution

It introduces a hybrid lattice Boltzmann simulation to analyze the effects of activity on highly asymmetric active emulsions, highlighting new morphological transitions and defect dynamics.

Findings

Activity reduces defects in extensile emulsions at low levels.

High activity induces aster-like rotating droplets in extensile systems.

Contractile activity leads to elongated structures at high levels.

Abstract

We investigate numerically, by a hybrid lattice Boltzmann method, the morphology and the dynamics of an emulsion made of a polar active gel, contractile or extensile, and an isotropic passive fluid. We focus on the case of a highly off-symmetric ratio between the active and passive components. In absence of any activity we observe an hexatic-ordered droplets phase, with some defects in the layout. We study how the morphology of the system is affected by activity both in the contractile and extensile case. In the extensile case a small amount of activity favors the elimination of defects in the array of droplets, while at higher activities, first aster-like rotating droplets appear, and then a disordered pattern occurs. In the contractile case, at sufficiently high values of activity, elongated structures are formed. Energy and enstrophy behavior mark the transitions between the different regimes.

Tables

Table 1. Table 1: Equilibrium distribution coefficients

$A_0=\rho -20A_2$	$A_1=4A_2$	$A_2=\frac{\left[{\Pi }_{\alpha \beta }+𝐯\zeta \left(Q_{\alpha \beta }+\frac{1}{3}{\delta }_{\alpha \beta }\right)\right]{\delta }_{\alpha \beta }}{24{𝐯}^2}$
$B_0=0$	$B_1=4B_2$	$B_2=\rho /12{𝐯}^2$
$C_0=-2\rho /3{𝐯}^2$	$C_1=4C_2$	$C_2=-\rho /24{𝐯}^2$
$D_0=0$	$D_1=8D_2$	$D_2=\rho /8{𝐯}^2$
$G_{0\alpha \beta }=0$	$G_{1\alpha \beta }=4G_{2\alpha \beta }$	$G_{2\alpha \beta }=\frac{\left[{\Pi }_{\alpha \beta }-\frac{1}{2}{\Pi }_{\delta \delta }{\delta }_{\alpha \beta }+𝐯\zeta \left(Q_{\alpha \beta }+\frac{1}{3}{\delta }_{\alpha \beta }\right)\right]}{8{𝐯}^4}$

Equations

$$\rho \left(\frac{\partial }{\partial t}+\mathbf{v}\cdot \nabla \right)\mathbf{v}$$

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot \left(\varphi \mathbf{v}\right)$$

$$\frac{\partial \mathbf{P}}{\partial t}+\left(\mathbf{v}\cdot \nabla \right)\mathbf{P}$$

$${\underline{\underline{\sigma }}}^{\text{total}}={\underline{\underline{\sigma }}}^{\text{passive}}+{\underline{\underline{\sigma }}}^{\text{active}}\text{mbox}.$$

$${\underline{\underline{\sigma }}}^{\text{passive}}={\underline{\underline{\sigma }}}^{\text{viscous}}+{\underline{\underline{\sigma }}}^{\text{elastic}}+{\underline{\underline{\sigma }}}^{\text{interface}}\text{mbox}.$$

$${\sigma }_{\alpha \beta }^{elastic}=$$

$${\sigma }_{\alpha \beta }^{interface}=\left(f-\varphi \frac{\delta F}{\delta \varphi }\right){\delta }_{\alpha \beta }-\frac{\partial f}{\partial \left({\partial }_{\beta }\varphi \right)}{\partial }_{\alpha }\varphi .$$

$${\sigma }_{\alpha \beta }^{active}=-\zeta \varphi \left(P_{\alpha }{P}_{\beta }-\frac{1}{3}|\mathbf{P}{|}^2{\delta }_{\alpha \beta }\right),$$

$$F$$

$$\psi_{6i}=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\exp{6\mbox{i}\theta_{ij}}\mbox{,}\qquad\Psi_{6}=\frac{1}{N}\biggl{|}\sum_{i=1}^{N}\psi_{6i}\biggr{|}\mbox{,}$$

$$\begin{array}{rl}& f_i(\mathbf{r}+{\mathbf{e}}_i\Delta t,t+\Delta t)-f_i(\mathbf{r},t)\\ & =-\frac{\Delta t}{\tau }[f_i(\mathbf{r},t)-f_i^{\text{eq}}(\mathbf{r},t)],\end{array}$$

$$\begin{array}{rl}& \sum _i{f}_i^{\text{eq}}=\rho \\ & \sum _i{f}_i^{\text{eq}}{e}_{i\alpha }=\rho {\mathbf{v}}_{\alpha }.\end{array}$$

$$\sum _i{f}_i^{\text{eq}}{e}_{i\alpha }{e}_{i\beta }={\Pi }_{\alpha \beta }+\rho {\mathbf{v}}_{\alpha }{\mathbf{v}}_{\beta },$$

$$\begin{array}{rl}{\Pi }_{\alpha \beta }& =-p{\delta }_{\alpha \beta }+2\xi (Q_{\alpha \beta }+\frac{1}{3}{\delta }_{\alpha \beta })Q_{\gamma ϵ}{H}_{\gamma ϵ}\\ & -\xi H_{\alpha \gamma }(Q_{\gamma \beta }+\frac{1}{3}{\delta }_{\gamma \beta })-\xi (Q_{\alpha \gamma }+\frac{1}{3}{\delta }_{\alpha \gamma })H_{\gamma \beta }\\ & -{\partial }_{\alpha }{Q}_{\gamma \nu }\frac{\delta f}{\delta ({\partial }_{\beta }{Q}_{\gamma \nu })}+Q_{\alpha \gamma }{H}_{\gamma \beta }-H_{\alpha \gamma }{Q}_{\gamma \beta },\end{array}$$

$$\begin{array}{rl}{f}_k^{\text{eq}}=& A_k+B_k{\mathbf{v}}_{\alpha }{e}_{k\alpha }+C_k{\mathbf{v}}^2+\\ & D_k{\mathbf{v}}_{\alpha }{\mathbf{v}}_{\beta }{e}_{k\alpha }{e}_{k\beta }+G_{k\alpha \beta }{e}_{k\alpha }{e}_{k\beta },\end{array}$$

Full text

1 Introduction

The capability of different systems of using energy taken from their environment to go out of thermal equilibrium, gives rise to a wealth of behaviors [1]. They range from swarming, self-assembly, spontaneous flows to other collective properties [2, 3, 4, 5]. This boosted a deep interest in addressing their study in order to look for possible new physics, explore common features between different systems, and develop new strategies in designing synthetic devices and materials with smart properties.

Self-propelled objects represent a remarkable example of active matter. Starting from the seminal model of Vicsek [6] for swarms, it was later realized that common features can be traced in several systems at different scales promoting the introduction of statistical models able to describe such behaviors [7, 8, 9, 10]. Another example of active matter, sharing many properties with suspensions of swimmers, is made by active gels that have been introduced to model mixtures of polar biological filaments with (active) motor proteins [11, 12, 13]. Their continuum modeling is based on the liquid-cristal description of long filaments in the nematic phase supplemented with additional contributions to introduce the motor activity [14, 1, 2].

Research in this field has been mainly focused on single-component active systems and to a lesser extent on the behavior of solutions of active and passive components. Mixtures of self propelled and passive particles have been studied by Brownian-like simulations [15, 16], focusing on the role of activity in separating the two components of the mixtures. Binary fluids with an active component have been considered in [17, 18] showing that the active part may cause instabilities on active passive interface. Recently [19] a model has been introduced where emulsification of the active component is favored by the presence of surfactant added to the mixture. This model generalizes the aforementioned active gel theory to describe the behavior of a mixture of isotropic passive and polar active fluids. The goal was to have a system with a tunable amount of active material that can be dispersed homogeneously in the fluid. This would also represent a further important step in the study of active turbulent fluids [20, 21] (with the possibility of tuning the intensity but also the spatial distribution of energy input in the system).

In [19] a symmetric mixture of active and passive components was considered. Their equilibrium configurations are dominated by the formation of local ordered lamellae [22, 23, 24, 25, 26]. It was shown that that activity may modify such a configuration leading to a variety of morphologies whose formation strongly depends on the intensity and the kind of active doping. Indeed, polar active fluids are said to be either extensile (e.g. bacterial colonies and microtubules bundles) or contractile (e.g. actine and myosin filaments) according to the nature of the stress exerted by the active component on its neighborhood. Furthermore, intensity of active doping can be tuned by keeping under control the amount of fuel available to active particles. This corresponds experimentally to keeping under control the amount of ATP in active gels of bundled microtubles [27] or the amount of oxygen available, the concentration of ingredients, or the temperature in bacterial suspensions. In the present model this is done by introducing a parameter representing the strength of the active stress acting in the system (see Section 2). The main result was that, even if under symmetric conditions, activity modifies lamellar configurations into an emulsion of passive droplets in an active matrix at sufficiently high contractile activity. On the other hand, lamellae change their morphology into rotating active droplets in the extensile case.

In this work we complement the previous analysis by considering a highly off-symmetric mixture with a $10:90$ ratio of the active and passive components both for extensile and contractile systems. Here the equilibrium state of the fluid is characterized by an ordered array of droplets of the minority phase positioned at the vertices of a triangular lattice. We will show that, despite the strong unbalance between the two components, activity greatly affects the morphology of the system, leading to the development of a wide range of patterns both for the concentration and the velocity field. In the extensile case a small amount of activity favors the elimination of defects in the system, as shown by measuring the number of defects in Voronoi tessellation. By increasing activity, isolated droplets tend to merge forming larger ’islands’ of active material and then, at still larger activity, big rotating droplets are observed. In the contractile case activity promotes the rupture of the hexagonal phase and the appearance of a matrix of the active component in the passive flowing background, differently from what happened in the symmetric case. The morphological study is supported by the analysis of the kinetic energy and enstrophy behavior.

The dynamic equations for the concentration of the active material and the polarization, which fixes the average orientation of the active component, are derived from a proper free-energy functional and supplemented with the Navier-Stokes equations for the whole fluid in the incompressible limit. These equations are numerically solved by using a hybrid lattice Boltzmann method [17, 28, 29, 30, 31, 32].

The paper is organized in the following way. The next Section is devoted to present the thermo-hydrodynamic description of the system, the numerical method, and the set of parameters relevant for the present study. In Sections 3 and 4 the results for the morphology and the corresponding flow patterns will be shown and related to the observed behaviors of energy and enstrophy. Finally, a discussion with some remarks and possible future lines of investigation will conclude the paper.

2 Model

We outline here the hydrodynamic model and the numerical method used to conduct our study. We consider a fluid comprising a mixture of active material and solvent with total mass density $\rho$. The physics of the resulting composite material can be described by using an extended version of the well-established active gel theory [1, 2, 17, 32, 33, 34]. The hydrodynamic variables are the density of the fluid $\rho$, its velocity ${\bf v}$, the concentration of the active material $\phi$, and the polarization ${\bf P}$, which determines the average orientation of the active material. The dynamic equations ruling the evolution of the system are

[TABLE]

in the limit of incompressible fluid. The first one is the Navier-Stokes equation, where $p$ is the isotropic pressure and $\underline{\underline{\sigma}}^{total}$ is the total stress tensor [35]. Eqs. (2)-(3) govern the time evolution of the concentration of the active material and of the polarization field, respectively. Since the amount of active material is locally conserved, the time evolution of the concentration field can be written as a convection-diffusion equation, Eq. (2), where $M$ is the mobility, $F$ a free energy functional as defined later, and $\delta F/\delta\phi$ is the chemical potential. The dynamics of the polarization field follows an advection-relaxation equation, Eq. (3), borrowed from polar liquid crystal theory. Here $\Gamma$ is the rotational viscosity, $\xi$ is a constant controlling the aspect ratio of active particles (positive for rod-like particles and negative for disk-like ones), $h=\delta F/\delta\textbf{P}$ is the molecular field. $\underline{\underline{D}}=(\underline{\underline{W}}+\underline{\underline{W}}^{T})/2$ and $\underline{\underline{\Omega}}=(\underline{\underline{W}}-\underline{\underline{W}}^{T})/2$ represent the symmetric and the antisymmetric parts of the velocity gradient tensor $W_{\alpha\beta}=\partial_{\beta}v_{\alpha}$, where Greek indexes denote Cartesian components. These contributions are in addition to the material derivative as the liquid crystal can be rotated or aligned by the fluid [35]. The stress tensor $\underline{\underline{\sigma}}^{total}$ considered in the Navier-Stokes equation of the model, Eq. (1), is splitted into the equilibrium/passive and non-equilibrium/active part:

[TABLE]

The passive part represents elastic response from solutes and is, in turn, the sum of three terms:

[TABLE]

The first term is the viscous stress written as $\sigma_{\alpha\beta}^{viscous}=\eta(\partial_{\alpha}v_{\beta}+\partial_{\beta}v_{\alpha})$ where $\eta$ is the shear viscosity. The second term is the elastic stress analogous to the one used in nematic liquid crystal hydrodynamics  [35]:

[TABLE]

where $\kappa$ is the elastic constant of the liquid crystal and the parameter $\xi$ depends on the geometry, as already mentioned. In addition $\xi$ establishes whether the fluid is flow aligning ($|\xi|>1$) or flow tumbling ($|\xi|<1$) under shear. The third term is borrowed from binary mixtures theory. It includes interfacial contribution between the passive and the active phase:

[TABLE]

Here $f$ is the free energy density. The active contribution to the stress tensor, not stemming from the free energy, is given by [33, 36]

[TABLE]

where $\zeta$ is the activity strength that is positive for extensile systems (pushers) and negative for contractile ones (pullers). The active stress drives the system out of equilibrium by injecting energy into it and satisfies the symmetry ${\bf{P}}\rightarrow-{\bf P}$.

The thermodynamics properties of the binary mixture, in absence of activity, are encoded in the following free-energy functional that couples the Landau-Brazovskii model  [37] to the distortion free-energy of a polar system:

[TABLE]

This is a generalization of the free energy functional for active binary mixtures defined in [33]. The first term, multiplied by the phenomenological constant $a>0$, describes the bulk properties of the fluid, the second and third ones determine the interfacial tension. Notice that here a negative value of $k$ favors the formation of interfaces while a positive value of $c$ has to guarantee the stability of the free-energy [37]. The Landau-Brazovskii model, with only the $\phi$ terms in the first line of Eq. (9), when the composition is symmetric, has a transition line for $k_{cr}=-\sqrt{8c\tilde{a}/15}$, in the mean field approximation [38], where $\tilde{a}=a\phi_{0}^{2}/8\phi_{cr}^{4}$ is half the coefficient of the $\phi$ quadratic term in (9). Lowering $k$ from positive to negative values leads the system to move from pure ferromagnetic phase to configurations where interfaces between components are favored; for values lower than $k_{cr}$ the system exhibits a periodic behavior such that the equilibrium state of the system is characterized by lamellae. For asymmetric compositions droplets of the minority phase are stable [22, 39, 40]. The bulk term is chosen in order to create two free energy minima, one ($\phi=0$) corresponding to the passive material and the other one ($\phi\simeq\phi_{0}$) corresponding to the active phase; $\phi_{cr}=\phi_{0}/2$, where $\phi_{cr}$ is the critical concentration for the transition from isotropic ($|{\bf P}|=0$) to polar ($|{\bf P}|>0$) states.

The bulk properties of the polar liquid crystal are instead controlled by the $|\textbf{P}|^{2}$ and $|\mathbf{P}|^{4}$ terms, multiplied by the positive constant $\alpha$. The choice of $\phi_{cr}$ has been made in order to break the symmetry between the two phases and confine the polarization field in the active phase $\phi>\phi_{cr}$. The second last term proportional to $(\nabla\mathbf{P})^{2}$ describes the energy cost due to elastic deformation in the liquid crystalline phase, gauged by the elastic constant $\kappa$ (in the single elastic constant approximation). Finally, the last term takes into account the orientation of the polarization at the interface of the fluid. If $\beta\neq 0$, P preferentially points perpendicularly to the interface (normal anchoring): towards the passive (active) phase if $\beta>0$ ($\beta<0$).

The equations of motion of the exotic polar active emulsion, Eqs. (1)-(3), are solved by means of a hybrid lattice Boltzmann (LB) scheme, which combines a LB treatment for the Navier–Stokes equation (see A for more details) with a finite-difference predictor-corrector algorithm to solve the order parameter dynamics.

Simulations have been performed on a periodic square lattice of size $L=256$. The concentration $\phi$ ranges from $\phi\simeq 0$ (passive phase) to $\phi\simeq 2$ (active phase). Unless otherwise stated, parameter values are $a=4\times 10^{-3}$, $k=-6\times 10^{-3}$, $c=10^{-2}$, $\alpha=10^{-3}$, $\kappa=10^{-2}$, $\beta=0.01$, $\Gamma=1$, $\xi=1.1$, $\phi_{0}=2.0$, and $\eta=1.67$. All quantities in the text are reported in lattice units.

We have considered a $10:90$ mixture of the active and passive component following two initialization procedures. In one case we started from a random configuration of polarization and concentration. The latter has been initialized considering a random variation of ten percent around its average value. Starting from this configuration, the system has been equilibrated without activity. Then, activity has been switched on and the evolution of the system studied. In the other procedure the random initial configuration has been evolved in the presence of activity. We checked that the two procedures lead to the same behaviour at late times. The results presented in the following generally are obtained with the second initialization procedure.

3 Extensile case

3.1 Small activity, hexatic order, and defects

Having in mind to fully characterize the behavior of the system as the value of activity changes, we begin by considering the case when the activity is off. In this case, at equilibrium, the system is characterized by an ordered array of droplets as can be seen looking at the contour plot of the concentration $\phi$ in Fig. LABEL:img:no_activity_phi. The droplets (and their centers of mass) can be easily pinpointed by putting a proper cutoff on the concentration field to distinguish active regions from passive ones. Each closed region of lattice sites that falls beyond the cutoff, is identified as a droplet. A good choice for the cutoff with our choice of the parameters is seen to be $\phi\sim 1.5$. Droplets are hexatically ordered, that is they occupy vertices of a triangular lattice, besides the presence of some defects. A Voronoi tessellation is used in order to unambiguously identify the nearest-neighbors network for the centers of mass of each droplet. Voronoi tessellation establishes a partitioning of the space with one closed region for each center of mass, according to the following rule: the region associated to the $i$-th droplet contains all the points of the space that are closer to its center of mass than to any other droplets. In Fig. LABEL:img:no_activity_vor it is shown the use of this analysis. Droplets with 5 nearest neighbors are highlighted in red, while those with 7 neighbors in blue. Observe that most of the defects appear in pairs indicating the presence of dislocations in the system[41].

For non-zero, but still small (positive) values of $\zeta$, the hexatic-droplet phase survives. At sufficiently high values of $\zeta$, but still in the small activity regime, activity is able to reduce the number of defects, driving the system in an almost completely ordered phase, as shown in the time evolution of Fig. 2, at $\zeta=0.006$. The system first forms transient differently-oriented domains bounded by grain-boundary defects (an example of a closed grain boundary is shown in Fig. LABEL:img:vor_0_006_3), then they shrink during the time evolution and disappear at the end, leaving a single hole in the layout. We checked taht if instead of switching on activity from the beggining, we start from configurations equilibrated without activity, and then we switch it on, the defects dynamics is not significantly affected. In the course of this evolution we also observe other features, here only transient, that play a major role at larger activity. In Fig. LABEL:img:0_006_2 few larger aster-like droplets are observed. We refer as aster-like droplets, big rotating droplets or simply asters to non-circular droplets which have the shape of an aster associated to the formation of vortices in the velocity field, as will be seen later. The presence of these aster-like droplets will be a predominant feature in the large activity regime as will be seen further in the text. We observe also a slight increase of droplets size during time evolution of the system as highlighted in Fig. LABEL:img:0_006_3.

At the end, for $\zeta=0.006$ and slightly larger values of $\zeta$, activity is such to let the droplets rearrange in an almost defect-free final configuration. To better understand the effect of activity on the hexatic arrangement we plot in Fig. 3 the defects ratio $D$ (with respect to the number of droplets) as the activity varies. We see that for very small values of activity, the hexatic order survives but there is no appreciable change in the number of defects with respect to the passive limit. However, for values of $\zeta$ around $0.006$, or slightly larger, for the longest simulated times, the hexatic order is enhanced as activity is such to remove the defects. At larger values of $\zeta$ a change in the overall morphological behavior of the mixture will be observed and the defects analysis looses significance. In order to further characterize order in the droplet phase we use a hexatic order parameter, as the one used for studying hexatic order and phase transitions in two-dimensional systems [42]. Local hexatic order parameter and its global average are defined, as in [43], by

[TABLE]

where $N_{i}$ is the number of neighbors of the $i$-th droplet, $\theta_{ij}$ is the angle that droplets $i$ and $j$ form with a reference axis, and $N$ is the total number of droplets. As shown in the inset of Fig. 3, the measure of $\Psi_{6}$ is consistent with the picture found by looking at the number of defects. As activity is increased from [math], the global hexatic parameter fluctuates around a positive value, which means that the system is partially hexatically ordered except for some defects. Within the activity range where we see no defects, $\Psi_{6}$ reaches its saturation value 1, being the system in a perfectly ordered phase.

3.2 Larger activity and aster-like droplets

Starting from $\zeta=0.0072$ aster-like droplets are no more transient and the system morphology, at late times, becomes dominated by the presence of these structures111We checked the behavior until $t=3.5\times 10^{6}$ as can be seen from the time evolution of the kinetic energy in Fig.5. . The transition to this regime can be analyzed, qualitatively, looking at the panels in Fig. 4. Here a series of snapshots of $\phi$ contour plots for different values of activity are reported, starting from $\zeta=0.0072$ up to $\zeta=0.0090$. A small number of isolated aster-like droplets linger at late times for $\zeta=0.0072$ (Fig. LABEL:zeta=0.0072), while for bigger values of activity these structures start grouping (Fig. LABEL:zeta=0.0074) filling progressively a large portion of the system (Figs. LABEL:zeta=0.0078, LABEL:zeta=0.0080 and LABEL:zeta=0.0090) at the expense of small droplets.

The mechanism of formation of the aster-like droplets can be intuited looking at Figs. LABEL:zeta=0.0072 and LABEL:zeta=0.0074. First, bigger droplets thicken in a portion of the system, leading two of them to merge, then the new big droplet starts rotating while incorporating neighboring droplets until an aster-like structure is fully stable. The time evolution of the kinetic energy per unit volume $E_{k}=\sum_{\mathbf{r}}\rho\mathbf{v}(\mathbf{r})^{2}/2$, for $\zeta=0.0072$, is shown in Fig. 5 (blue curve) and reproduces the aster formation: The system sits in the droplet configuration with a constant value of kinetic energy density until the merging of the two droplets leads to a strong discontinuity in the energy density. When the aster-like droplets are fully stable, there is a strong contribution to the kinetic energy coming from their angular velocity. These structures then begin to slow down until they disappear, and the energy density returns to a lower constant value. Then, again, new asters are formed and do not disappear even at late times. In contrast, at lower activity (yellow curve in Fig. 5, for $\zeta=0.006$) asters disappear and the kinetic energy goes to zero remaining constant for the time checked. The red curve in Fig. 5 corresponds to the case $\zeta=0.0078$. The kinetic energy is appreciably greater then the previous cases and increases until a stationary state is reached. For this value of activity asters are continuously formed and never disappear in the system.

The velocity field at different activities is displayed in Fig. 6. We see that the resulting aster-like droplets (Fig. LABEL:img:zoom0.0072) have a velocity field much stronger that the one around droplets (see caption of Fig. LABEL:img:zoom0.0072). They behave as a vortical source for the velocity field and for this reason we refer to them also as big rotating droplets. The velocity field as activity varies reproduces the behaviors so far discussed. At $\zeta=0.006$ (Fig. LABEL:img:zoom0.006) we distinguish vortical structures that are not localized anymore in correspondence of the droplets, as for smaller values of the activity (Fig. LABEL:img:zoom0.0030), but start mixing with neighboring vortices.

3.3 Very large activity and overview

At even larger values of activity ($\zeta\gtrsim 0.01$) the aster-like structures are destroyed by the flow and no clear morphological pattern can be observed. The system appears completely mixed as it will be shown.

A better characterization of the behavior of the system in this regime comes from the study of hydrodynamic quantities such as the kinetic energy per unit volume $E_{k}$ and the enstrophy per unit volume defined as $\varepsilon=\sum_{\mathbf{r}}\mathbf{\omega}(\mathbf{r})^{2}/2$, where $\omega_{i}=\epsilon_{ijk}\partial_{j}\mathrm{v}_{k}$ is the vorticity vector and $\epsilon_{ijk}$ is the completely anti-symmetric Levi-Civita tensor. Kinetic energy represents a measure of the strength of fluid flows in the system, while enstrophy can be used to figure out whether the velocity field has developed vortical structures.

The graphs in Fig. 7 show the behavior of energy and enstrophy averaged in time over uncorrelated configurations as functions of the activity. They are both null for $\zeta=0$ and raise to reach a peak in the range $0.01-0.02$. Values are significantly different from [math] at $\zeta=0.0072$ and this corresponds to the formation of vortical asters. If activity is raised over $\zeta\simeq 0.007$, the number of vortical structures increases until $\zeta\simeq 0.02$. For larger values of $\zeta$ asters become unstable and flexible, elongated rotating structures are formed until one arrives at configurations like that of Fig. 8, at $\zeta=0.09$, with the velocity field exhibiting a pattern like the one shown. When these structures start melting, rotational contribution to the kinetic energy starts decreasing while the kinetic energy stays approximately constant. In this regime weak vortical flows span the system. Such flows dissipate energy while moving, according to the dissipative nature of the fluid, and weaken in intensity until they expire in small and slow vortices or simply merge in more intense flows as suggested by the streamtraces displayed in Fig. 8.

4 Contractile case

We now consider the case of contractile ($\zeta<0$) emulsions. For slightly negative values of $\zeta$, configurations at late times are similar to the extensile case.

Increasing $\lvert\zeta\rvert$, the hexatic order is progressively lost: First a change in the droplet dimension and shape is observed, as it can be seen in Fig. LABEL:img:-0.003. This leads the fluid flow to assume different configurations throughout the system. At very low contractile activity, not shown, the velocity field is small in magnitude and enforced in regions close to droplets creating vortical structures that rotate around the center of the droplets themselves. For values of activity, such as $\zeta=-0.008$, droplets start oscillating, following the flow (Fig. LABEL:img:zoom-0.003) and, eventually, merging (see Figs .LABEL:img:-0.003 and LABEL:img:zoom-0.003) thus creating elongated structures. For large enough activity the hexagonal pattern is lost (Fig. LABEL:img:-0.013). Flow fields are sustained by the active stress that continuously drives energy into the system and generates flows that are not confined around the droplets anymore, but move throughout the system as it is shown in Fig. LABEL:img:zoom-0.02. This happens for values of $\zeta$ ranging from $-0.01$ to $-0.02$ but for even stronger values, the merging of structures throughout the system affects the equilibrium configuration so strongly that when $\zeta\lesssim-0.02$, it is not possible to distinguish any kind of definite pattern anymore.

Here again, together with the morphology it is useful to look at the behavior of the kinetic energy and enstrophy per unit volume. The graphs in Fig. 10 show their behaviors averaged in time over uncorrelated configurations of the system for different values of the activity $\zeta$. Both energy and enstrophy are null when $\zeta=0$, according to the fact that in absence of an energy source no fluid flow can be sustained in a dissipative fluid; then they increase until $\zeta\simeq-0.01$ when enstrophy reaches a peak. This trend fits the morphology behavior presented before: As far as droplets are preserved, the velocity field creates vortical structures around them that raise in strength for stronger values of activity, but when the merging of droplets starts affecting the morphology, the local vortical structures are progressively replaced by flows without ordered structure. This brings energy to stay approximatively constant for values $\zeta\lesssim-0.02$ while enstrophy collapses rapidly.

5 Conclusions and perspectives

This work shows how activity modifies the equilibrium hexatic-ordered droplets configuration of an highly asymmetric active emulsion.

In the extensile case we found three different regimes. For $\zeta$ up to $0.005$ hexatic order is preserved, while is enhanced for values of $\zeta$ around $0.006$, which is the critical value for a transition to a regime dominated by big-rotating droplets. In fact, starting from $\zeta=0.0072$ the system morphology becomes characterized by the presence of aster-like droplets. At even larger values of activity ($\zeta\gtrsim 0.01$) the aster-like structures are destroyed by the flow and the system appears completely mixed. We have scored the transitions studying energy and enstrophy behaviors. In the case of contractile emulsions, for slightly negative values of $\zeta$ stationary configurations are similar to the extensile case while, strong contractile activity gives rise to elongated structures.

In the future we plan to extend these studies to three dimensional systems, where even richer morphologies ad flow patterns ca be expected.

6 Acknowledgments

Simulations have been performed at Bari ReCaS e-Infrastructure funded by MIUR through the program PON Research and Competitiveness 2007-2013 Call 254 Action I. We thank D. Marenduzzo, A. Tiribocchi and F. Bonelli for the very useful discussions.

Appendix A Numerical method

In the LB scheme the evolution of the fluid is defined in terms of a set of $N$ discrete distribution functions $\{f_{i}\}$ $(i=0,...,N-1)$ which obey the dimensionless Boltzmann equation in the BGK approximation:

[TABLE]

where $\mathbf{r}$ and $t$ are the spatial coordinates and the time, respectively, $\{\bf{e_{i}}\}$ $(i=0,...,N-1)$ is the set of discrete velocities, $\Delta t$ is the time step, and $\tau$ is a relaxation time which characterizes the relaxation towards the equilibrium distributions $f_{i}^{\textit{eq}}$. The shear viscosity $\eta$ is related to $\tau$ by the relationship $\eta=c^{2}\Delta t\frac{2\tau-1}{6}$. The value of $N$ depends on the space dimensions and the lattice geometry. The moments of the distribution functions allow to write the conservation laws for the density and total momentum in the form:

[TABLE]

The second moment, which describes the balance between energetic densities at stake, is fixed in order to find the hydrodynamic equations in the continuum limit. It is given by

[TABLE]

where $\Pi_{\alpha\beta}$ represents the pressure tensor. Introducing the nematic tensor $Q_{\alpha\beta}={\bf P}_{\alpha}{\bf P}_{\beta}-\frac{1}{3}|{\bf P}^{2}|\delta_{\alpha\beta}$, the pressure tensor can be written as

[TABLE]

where $H$ is the field conjugated to the nematic tensor and $p$ the ideal gas pressure.

The equilibrium distributions are expanded up to the second order in the velocities

[TABLE]

where the index $k$ labels the different directions on the discretized lattice. The $A_{k}$, $B_{k}$, $C_{k}$, $D_{k}$ and $G_{k\alpha\beta}$ are characteristic parameters to be determined to get the right hydrodynamic equations in the continuum limit by imposing the conditions given by Eqs. (12)-(13). For a two-dimensional square lattice with $N=9$ velocities (D$2$Q$9$), which is the model here considered, these parameters are reported in Table 1, and the lattice velocities are ${\bf e}_{0}=(0,0)$, ${\bf e}_{1,2}=(\pm u,0)$, ${\bf e}_{3,4}=(0,\pm u)$, ${\bf e}_{5,6}=(\pm u,\pm u)$, ${\bf e}_{7,8}=(\mp u,\pm u)$ with $u=\Delta x/\Delta t$, $\Delta x$ being the lattice mesh size. In Table 1 $k=0$ corresponds to the rest lattice velocity, while $k=1$ and $k=2$ correspond to the velocities directed towards the first and second neighbors, respectively.

The LB scheme is coupled with a finite-difference predictor-corrector algorithm, in order to solve simultaneously Eqs. (1)-(3). This method was previously used to study the hydrodynamics of binary fluids [28, 44], liquid crystals and active matter [30, 45].