Shearing effects on the phase coarsening of binary mixtures using the Active Model B

The phase separation of a two-dimensional active binary mixture is studied under the action of an applied shear through numerical simulations. It is highlighted how the strength of the external flow modifies the initial shape of growing domains. The activity is responsible for the formation of isolated droplets which affect both the coarsening dynamics and the morphology of the system. The characteristic dimensions of domains along the flow and the shear direction are modulated in time by oscillations whose amplitudes are reduced when the activity increases. This induces a broadening of the distribution functions of domain lengths with respect to the passive case due to the presence of dispersed droplets of different sizes.

The phase separation of a two-dimensional active binary mixture is studied under the action of an applied shear through numerical simulations. It is highlighted how the strength of the external flow modifies the initial shape of growing domains. The activity is responsible for the formation of isolated droplets which affect both the coarsening dynamics and the morphology of the system. The characteristic dimensions of domains along the flow and the shear direction are modulated in time by oscillations whose amplitudes are reduced when the activity increases. This induces a broadening of the distribution functions of domain lengths with respect to the passive case due to the presence of dispersed droplets of different sizes.

I. INTRODUCTION
Over the last decades much research has been dedicated to studying active materials, i.e., non-equilibrium systems in which internal units continuously consume energy, usually stored in the environment, to self-propel [1][2][3][4][5]. Examples of active matter range from bacterial and algal suspensions [6][7][8][9] to colloidal particles acquiring motion through chemical reactions [10][11][12], up to the cytoskeleton of living cells [13][14][15], to name but a few. Importantly, to convert their energy into motion, particles violate time-reversal symmetry (TRS) at the micro-scale.
An active process of particular relevance to us is the motility-induced phase separation (MIPS), where a suspension of repulsive motile particles can phase separate into bulk dense (liquid) and dilute (vapor) regions [16][17][18]. This essentially occurs because active particles aggregate where they move slowly, thus leading to a local density increase causing a slowing down and further accumulation. The resulting clusters show highly dynamical particle exchange leading to macroscopic phase separation [19][20][21]. The kinetics of MIPS has been found to share many features with that of passive systems with attractive interactions, thus suggesting that TRS could be restored at the macroscale level [16,17,[22][23][24][25][26][27]. However, large scale simulations of active Brownian particles [20,21,28,29] as well as experiments [10][11][12] reveal that MIPS may also exhibit further non-equilibrium features, such as the formation of mesoscopic vapor bubbles within the particles aggregates and microphase separation [20,[28][29][30]. This scenario would indeed support the view that time-reversal symmetry is manifestly broken macroscopically.
Alongside particle-based simulations, the phenomenology of MIPS can be also described using continuum field theories, either by an explicit coarse-graining of the microscopic dynamics or via symmetry arguments and conservation laws [31]. The latter strategy has led to the construction of the Active Model B (AMB) [24,28,32], a field theory that, in the absence of solvent, essentially extends the passive model B [33] by incorporating an active gradient term that cannot be obtained from a free-energy functional, thus breaking time-reversal symmetry. The passive Model B (MB) falls within a class of phase-field models describing the kinetics of phase separation through a conserved scalar order parameter φ capturing, for example, the density of colloids or the concentration in a binary mixture. If hydrodynamics can be neglected, the evolution of φ is governed by a Cahn-Hilliard equation [34,35], where the thermodynamic force driving the relaxation of φ stems from the functional derivative of a φ 4 free energy with square gradient terms. Such simplified description captures, for example, the result R ∼ t α , for the dependence of the domain size R on time t with a characteristic growth exponent α = 1/3 [36].
MB provides a robust framework to describe the passive phase separation subject to an external shearing [37][38][39][40][41][42][43]. In that case the growth is anisotropic, being characterized by domains stretched and tilted along the flow direction as observed in simulations [42,44] and experiments [45]. Indeed, numerical solutions of the dynamical equation with shear [42] confirm the presence of anisotropy of growing domains, which are found to exhibit two typical lengths. Theoretical calculations based on a renormalization group analysis shows that the growth exponent α y in the shear direction, perpendicular to the flow, is not changed by the applied velocity profile, while the exponent α x along the flow direction is increased by 1 [41]. Moreover, relevant physical quantities are found to oscillate on a logarithmic time scale, due to a periodic stretching and breaking-up of domains. The interest towards the role of shear in phase ordering is still vivid, as witnessed by recent studies [46,47].
In the context of active matter, while efforts have been addressed to investigate the active phase separation in the absence of an external driving [20-22, 28, 48, 49], much less is known about the dynamic response observed when a shear flow is applied. Though it was shown [50] that in AMB the transition from homogeneous to bulk phase separation belongs to the same universality class of equilibrium MB, numerical studies of coarsening in AMB provide different scenarios. In Refs. [32,51], for example, it was argued that in AMB without shear the size of growing domains is still compatible with a growth exponent α = 1/3 (i.e., activity has negligible effects on the coarsening dynamics) with a change in the static phase diagram owing to a pressure jump, even in the case of a flat interface. On the contrary, very recent simulations [52,53] put forward the existence of a late-time growth exponent α = 1/4, a result akin to that obtained using particle-based simulations [20,54].
In this work we aim at characterizing, in the AMB theory, the shear-induced morphology of growing domains in two spatial dimensions focusing on the role played by the external flow and the activity. Our results show that, in agreement with the AMB without shear, the activity promotes the formation of isolated drops. This significantly affects the dynamic behavior of the mixture subject to a shear flow. Indeed, unlike the passive counterpart, the amplitude of the time oscillations of domain sizes diminishes for increasing values of activity, although their growth exponents are only weakly modified. In addition, the probability distribution functions of the lengths of patterns computed along the two spatial directions are found to broaden as the activity augments.
The paper is structured as follows. In Section II the model is properly defined. The numerical results are presented and discussed in Section III. The phenomenology of phase separation is discussed in weak and strong regimes, remarking the differences with respect to the passive model. The role played by isolated droplets appearing during coarsening is highlighted in connection with the time evolution of the size of domains along the two directions. Finally, some conclusions are drawn.

II. THE MODEL
We consider a system with a scalar order parameter φ(r, t) at position r and time t. In the framework of the present model, the order parameter can be considered as the deviation of the density with respect to a reference value. It satisfies a conserved dynamics given by the following equation Thermal fluctuations are not considered here as it is usually done when studying phase separation [36]. At the l.h.s. a convective term couples φ to an imposed linear shear flow. This has the form v =γye x whereγ is the shear rate, y is the coordinate along the (shear) y−direction and e x is the unit vector along the (flow) x−direction. The addition of this term to the original equation of the AMB is the main novelty of this study which allows us to consider the effect of an external velocity profile on the active coarsening dynamics. Hydrodynamic effects are not taken into account since the evolution Equation (1) is not coupled to the Navier-Stokes equation. This is because we are interested in considering only diffusive phase separation. Indeed, to account for the fluid motion in an active medium, the dynamics of φ needs to be coupled to that of a momentum-conserving solvent, whose evolution is governed by the Navier-Stokes equation including an active contribution in the stress tensor [55].
The current J is proportional to the negative gradient of a chemical potential where the mobility Γ is set to unity in the following. The chemical potential µ = µ P + µ A is the sum of passive and active contributions given, respectively, by The term (3), where a, b, κ are positive and of the order unity, corresponds to the passive contribution of the model B, and can be obtained from the derivative of the square gradient ϕ 4 free-energy functional. In the passive case, the system would demix in two coexisting states at binodal densities φ P 1,P 2 = ±φ P with φ P = a/b, corresponding to the minima of the free-energy density. The energy cost for the formation of interfaces is proportional to κ. The active contribution (4), on the contrary, cannot be obtained from a proper free energy and is the simplest choice at second order in gradients. It was introduced in the so-called active model B [32] in order to explicitly break TRS. For completeness we note that this term is similar to the one appearing in other nonlinear partial differential equations. Among others we cite, e.g., the Kardar-Parisi-Zhang equation for nonlinear interfacial diffusion [56], the Hunter-Saxton equation for nematic liquid crystals [57], and the Kuramoto-Sivashinsky equation for instabilities in laminar flame front [58,59]. Such a term is controlled by the active parameter λ which can be adjusted to go from the passive case (λ = 0) to the active one (λ = 0) with λ ∼ O(1) [32]. We remark that the presence of the λ term is such that the system is invariant under the transformation (φ, λ) → −(φ, λ). For this reason we restrict our attention to the case with λ ≥ 0. For the active model B the binodals φ 1,2 (φ 1 < 0 < φ 2 ) depend on λ and can be calculated by using the method put forward in Ref. [32] and generalized in Ref. [51]. The spinodals do not depend on λ [51] and are located at φ S1,S2 = ±φ S with φ S = a/(3b) as in the model B. The state with uniform density is locally unstable between the spinodals while it is metastable between the spinodals and the binodals. When λ increases, the binodal gets closer to the spinodal on the negative φ side and stays far on the other side.
Equation (1) is solved in two dimensions by using a finite-difference scheme. The field φ is discretized on the nodes (x i , y j ) (i, j = 1, 2, ..., L) of a square lattice with L × L nodes and mesh size ∆x. The time is discretized in time steps ∆t with time values given by t n = n∆t, n = 1, 2, 3, .... Any discretized function h at time t n on a node (x i , y j ) (i, j = 1, 2, ..., L) of the lattice is denoted by h(x i , y j , t n ) = h n ij . At each time step we update φ n → φ n+1 using an explicit first-order Euler algorithm for the time derivative [60] Central-difference schemes [61] are coded for the spatial derivatives. The x derivative is given by and analogously for the y derivative. The Laplacian operator appearing in the chemical potential µ P is discretized as Periodic boundary conditions (BC) are adopted along the flow direction and Lees-Edwards BC [62] are used in the shear direction. The latter take into account the space shiftγL∆x∆t occurring in a time step, due to shear, between the lower and upper row of the lattice. This is essentially done by identifying a point (x i , y 1 ) on the lower row of the lattice with the one placed on the upper row at (x i +γL∆x∆t, y L ) (i = 1, 2, ..., L).

III. RESULTS
In the following, all results are obtained with ∆x = 1 and ∆t = 0.001, values ensuring numerical stability. We checked that simulations are converged for these values of space and time discretization units by looking at the behavior of the binodal values and of averaged domain sizes for different initial conditions (see the following). We verified that results are stable upon decreasing the discretization units. To this purpose we considered some runs with ∆x = 0.5 and ∆t = 0.0005, and no significant differences were found in the binodals as well as in the typical extensions of forming patterns during coarsening. The model parameters in Equation (1) are a = b = 1/4 and κ = 1 while λ is varied within the range [0, 3]. The shear rateγ is changed to access different shear regimes in the phase-separation process. To this purpose we introduce a dimensionless shear rateγ =γt D where t D is the interface diffusion time. Weak and strong shear regimes are characterized by values ofγ lower and higher than unit, respectively [43].

A. Planar Interface
In the case of the model B one has t D = a is the width of the planar interface between two coexisting phases described by the function φ(x) = φ P tanh (2x/ξ), and σ = 2 3 φ 2 P √ 2κa is the interface tension [36]. For our choice of the parameters, it is t D = 384 in model units. However, for the active model B with λ > 0 there are no explicit expressions for ξ and σ as well as for the binodals, though it can be shown that there is a solution with a planar interface between the coexisting phases at densities φ 1,2 [32]. For this reason we first compute numerically the binodals φ 1,2 (λ) by considering the relaxation of a flat interface between two states initially set at the values ±φ P .
Once the values φ 1,2 (λ) are found, a sharp profile between the initial states set at φ 1,2 (λ) is let to evolve to the steady interface. The stationary profiles are shown in Figure 1 for the values λ = 0, 1, 2, 3. Numerical data are successfully fitted by a kink profile (see the Appendix for further details). The values of the interface width ξ obtained from fits for different values of the activity are presented in Figure 2 and show a linear dependence of ξ on λ for λ 1. The binodal densities φ 1 and φ 2 increase with λ as φ 1 approaches the spinodal φ S1 = −φ S , being φ S = √ 3/3 when a = b (like in our case). The diffusion time t D is computed as the time for the initial sharp profile to relax to the steady interface (see the Appendix IV). Numerical results give t D 4 × 10 2 at λ = 0 in agreement with the theoretical estimate of the model B. When λ > 0, it is found that t D increases with λ reaching the highest value t D 1.1 × 10 3 when λ = 3. The values of t D are plotted in the right panel of Figure 2 and show that t D increases linearly with λ up to λ 1.5, before slowing down its growth. In the following we consider a range of values of the shear rate betweeṅ γ w = 0.9 × 10 −3 andγ s = 3.6 × 10 −3 . The valuesγ w andγ s are such to access weak (γ w t D (λ) 1) and strong (γ s t D (λ) > 1) shear regimes, respectively, for all considered values of λ. This guarantees that the shear regime is not affected when changing activity while keeping fixed the strength of applied flow.

B. Phase Separation under Weak and Strong Shear
Now we move on to the study of the phase separation under an external shear flow by varying the activity parameter λ forγ =γ w ,γ s . We consider a system initially prepared in a symmetric disordered state with ϕ(r, 0) = ω where ω is a random number in the range [−0.01, 0.01]. This state corresponds to a critical composition of the system. The size of the lattice is L = 1024 and measures are averaged over five independent runs. The different values of the dimensionless shear rateγ influence the initial morphology of the forming domains. This can be seen in Figures 3  and 4 where snapshots of systems at consecutive times are shown for the cases at weak (γt D = 1.0) and strong shear (γt D = 4.0), respectively, with λ = 3. At weak shear, domains form and grow while the density φ attains the steady binodal values before shear effects come into play at values of strainγt > 1. This occurs since the relaxation time scale t D (λ) is smaller than the shear characteristic time 1/γ. Once isolated drops at density φ 1 are produced (γt = 1), the applied flow advects them along the x−direction and favours their merging, thus causing the formation of elongated domains along the y−direction (γt = 1.75). Afterwards the shear stretches these structures which are tilted as well, being characterized by different thicknesses. Isolated droplets at density φ 1 can be observed inside domains at density φ 2 (γt = 5.3). Later on, the shear further deforms and tilts the domains, which may eventually break up. After the overstretching, domains retract forming φ 2 phases with a larger thickness where isolated droplets of the other phase get trapped (γt = 19.4). Such stretching and bursting persist though this periodic behavior cannot be observed over very long periods of time due to the finite size of the system.
In the strong shear regime, the initial morphology is different with respect to the previous case, since the interface diffusion time is larger than the inverse of the shear rate. Indeed, while domains grow (as quantified later) and start to be deformed by the flow, the two phases are still far from the steady binodal densities (see the panel atγt = 2 in Figure 4), an effect promoting an initial thinning of domains along the shear direction. Afterwards, the elongations and ruptures of domains previously described take place once again, a phenomenon also observed in simulations of passive model B with strong shear [41,42]. However, unlike such cases, the main feature here is that, under shear, isolated droplets of the φ 1 phase survive and are dispersed in the φ 2 matrix (see in particular the last snapshots of Figures 3 and 4). This is not the case if λ = 0, as demonstrated in Figure 5 where instantaneous configurations of the passive case taken at equal times are shown.
A scenario akin to that just described in the AMB occurs when considering sheared binary mixtures with surface diffusion [64]. However, in the present work the physical explanation is different. The fraction β of the φ 1 phase is larger than the one of the φ 2 phase when λ > 0 despite the fact that we consider a critical system with a symmetric initial composition such that < φ >= 0 (the symbol < ... > denotes an average over the lattice). Since the dynamics described by Equation (1) is conserved, the average value of the order parameter φ is preserved and it has to be βφ 1 + (1 − β)φ 2 = 0. In a passive mixture the values of the binodals φ 1 , φ 2 are symmetric with respect to the value φ = 0 and it results β = 0.5. In the AMB the binodals are not symmetric anymore with respect to the value φ = 0 and one has β > 0.5, since |φ 1 | < φ 2 when λ > 0. As a result, the φ 1 phase is more abundant than the φ 2 one. In our study it is φ 1 −0.71 and φ 2 1.15 when λ = 3 so that β 0.62. Therefore, in AMB, it is the activity that produces an effective off-symmetric mixture though the initial state is not, in agreement with previous studies of AMB [32]. Such off-symmetric mixture consists of droplets of the majority (φ 1 ) phase dispersed in the minority (φ 2 ) one, a situation persisting under shear. We finally note that a different morphology would be observed in the coarsening of off-symmetric passive mixtures under shear with a similar ratio between the two phases [65]. In this case, at low shear rate small droplets evaporate due to Ostwald ripening, while for strong flows it is the minority phase to be dispersed in a large number of droplets.

C. Domain Size
The coarsening dynamics can be investigated by looking at evolution of the typical measures of domains. The sizes R x and R y along the flow and the shear directions, respectively, are shown in Figure 6 and are computed as the inverse of the first moments of the structure factor R x,y (t) = π dkC(k, t) dk|k x,y |C(k, t) . , t) is the structure factor averaged over different realizations of the system and φ(k, t) is the Fourier transform of the order parameter.    Table I. Fitted values of the amplitudes A for different values of the activity λ in the weak shear regime.
At weak shear, it results that R x R y untilγt 1, then R y grows initially faster than R x due to the shear-induced dragging of droplets along the flow direction. The typical size along the shear direction attains a maximum atγt 1.75 and then decreases, while the growth rate along the flow direction increases due to shear stretching. The quantity R y oscillates on a logarithmic time scale between minimum and maximum values corresponding to elongated and broken domains, respectively. The amplitude of such oscillations diminishes when increasing the activity parameter λ and is related to the existence of small droplets that cannot grow any further. The radii R x for different values of λ show a similar trend compatible with an asymptotic power-law growth R x ∼ A(λ)t αx . Fitting the power law for R x to numerical data shows that the amplitudes A decrease with the activity as illustrated in Table I, while the exponent α x 1 does not seem to be affected by λ.
Due to the limited size of the simulated system and the superimposed oscillations, it is difficult to estimate the growth exponent α y along the shear direction although, atγt 5, one observes a regime with α x − α y 1. This result is due to the advection of domains by the applied flow, and has been observed in the passive model B under shear as well [41,42].
At increasing values of shear rate, the growth along the flow direction proceeds faster than that in the shear direction, with a local maximum of R x atγt 3 corresponding to a minimum for R y , since domains are strongly deformed by the shear. Note that the height of the maximum of R x is reduced by increasing the activity, a feature associated to the presence of droplets which can be deformed only slightly by the flow. Whenγt 7, R x shows a power-law growth with an exponent α x 1.1 not depending on λ, while the size R y exhibits a periodic behavior on a logarithmic time scale with amplitude that shrinks with λ. By using power counting, it might be expected that the net effect of the convective term is to increase the growth exponent of the unsheared system by 1 in the flow direction. However, as previously discussed, in AMB with no external flow the value of the growth exponent has not been definitely determined yet. Our results seem to point towards a picture in which the activity has mild effect on the coarsening dynamics, as in the unsheared case [32], since we find that α x 1 as in the passive case and α x − α y 1. Additional simulations run on systems of size L = 512 do not provide further insights for evaluating α x and α y as well as for attempting a finite-size scaling analysis. Indeed, a more accurate estimate of the growth exponents would likely require much larger systems, thus dramatically increasing the computational resources necessary to simulate their late-time dynamics.
To elucidate the role of flow strength, simulations with shear rateγ = 2.8 × 10 −3 are also considered. This value is such thatγt D (λ) 1.1 for 0 ≤ λ ≤ 3 so that the system is in the strong regime. The time behavior of R x and R y is similar to the one shown in Figure 6, with no significant effect on the growth exponents. We only observe a reduction in the amplitudes of oscillations at the smallest values of activity along the x-direction, since domains are less affected by flow. Being less deformed, they can grow along the shear direction, a result witnessed by wider amplitudes of R y with respect to the case withγ = 3.6 × 10 −3 . These effects reduce when λ increases, since, as previously discussed, the dynamics is deeply modified by the presence of droplets for the highest value of activity.
The size distribution of domains can be analyzed by calculating the normalized probability distributions P (L x,y ), having domains of lengths L x and L y along the flow and the shear directions, respectively. By moving along all rows (the flow direction) of the computational domain, L x is computed as the size of unidimensional domains with equal composition (same sign of the density φ). From the registered values of L x , the function P (L x ) is derived. The same procedure is adopted along the columns (the shear direction) of the lattice to determine L y and then P (L y ). In Figures 7 and 8 we plot P (L x ) and P (L y ) for λ = 0, 3 in the weak and strong shear regime. At weak shear (see Figure 7), P (L x ) and P (L y ) show two peaks in the active case and a single one in the passive situation. The peaks at the smaller values of L x and L y (with L x L y 8) correspond to the presence of isolated droplets, while the other peaks indicate elongated domains. In particular, P (L x ) broadens when going from the maximum of R y atγt 1.75 to the minimum atγt 5.3 and, then, shrinks when R y grows again. In the passive counterpart, the distribution functions are less broad. When looking at the shear direction, in the active case the peak of P (L y ) at the larger value of L y oscillates in height, with a local maximum attained when R y is at the minimum (γt 5.3). Later P (L y ) broadens. In the passive case the distribution P (L y ) has a single peak (narrower than in the active case) that oscillates in height, following the cyclic dynamics of elongations and ruptures.   Moving to the strong shear regime (Figure 8), we observe that P (L x ) is characterized by the presence of two peaks both in the passive and in the active case, roughly at the same positions. When λ = 3, the peak at L x 10 is always the prevailing one, owing to the presence of droplets spanning the system. The value of L x , corresponding to the second peak, increases in time attaining its maximum atγt 17.4, when domains are highly stretched (R y is at its minimum). The subsequent breaking-up of domains determines the shift in position of the second peak to smaller values. Moreover, it can be seen that the distributions are broader in the active case. This feature also holds for   P (L y ), whose main peak oscillates in height during the time evolution. We note here that in the active system the first peak of P (L y ) at the smaller value of L y , corresponding to isolated droplets, is less pronounced with respect to that of the weak shear regime. This can be attributed to the fact that, as previously discussed, initially forming droplets are fused by shear before growing in size.

IV. CONCLUSIONS
To summarize, we have investigated the phase separation of an active binary mixture subject to an applied shear flow. For this purpose we numerically solved the phenomenological equation of the active model B [32], supplemented by a convective term that couples the order parameter to the external velocity field. The initial morphology depends on the strength of shear. Later, growing domains are elongated, tilted, and burst by the flow. This is reflected in the typical sizes of domains R x and R y along the two spatial dimensions, which appear to be modulated by oscillations on a logarithmic time scale. However, the presence of activity is such that the fraction of one phase is larger than the other one, despite the initial symmetric composition. This induces the presence of droplets of the more abundant phase which span the system and are responsible for the observed reduction of the amplitudes of the oscillations of R x and R y when increasing the activity parameter. Though the limited size of the simulated system does not allow an estimate of the growth exponents α x and α y along the flow and the shear directions, respectively, we find that α x − α y 1 as in the passive case [41,42]. The combined effect of activity and shear on the overall morphology has been studied by considering the probability distribution functions (PDFs) of the size of patterns along the two spatial directions. Our simulations suggest that such PDFs are characterized by a width that broadens with activity.
We hope that our results may stimulate further research on this system. It would be of interest, for example, understanding the role played by the shear in active binary mixtures where all terms breaking time-reversal symmetry to leading order in ∇ and φ are included [28], as well as how thermal noise is expected to impact on morphology and growth dynamics in the presence of an external driving.