Phase separation dynamics in deformable droplets

Phase separation can drive spatial organization of multicomponent mixtures. For instance in developing animal embryos, effective phase separation descriptions have been used to account for the spatial organization of different tissue types. Similarly, separation of different tissue types and the emergence of a polar organization is also observed in cell aggregates mimicking early embryonic axis formation. Here, we describe such aggregates as deformable two-phase fluid droplets, which are suspended in a fluid environment (third phase). Using hybrid finite-volume Lattice-Boltzmann simulations, we numerically explore the out-of-equilibrium routes that can lead to the polar equilibrium state of such a droplet (Janus droplet). We focus on the interplay between spinodal decomposition and advection with hydrodynamic flows driven by interface tensions, which we characterize by a Peclet number $Pe$. Consistent with previous work, for large $Pe$ the coarsening process is generally accelerated. However, for intermediate $Pe$ we observe long-lived, strongly elongated droplets, where both phases form an alternating stripe pattern. We show that these ``croissant'' states are close to mechanical equilibrium and coarsen only slowly through diffusive fluxes in an Ostwald-ripening-like process. Finally, we show that a surface tension asymmetry between both droplet phases leads to transient, rotationally symmetric states whose resolution leads to flows reminiscent of Marangoni flows. Our work highlights the importance of advection for the phase separation process in finite, deformable systems.


Introduction
2][3] During this process, phase domains of different composition first emerge either through spinodal decomposition or nucleation, before these domains coarsen through diffusive processes such as Ostwald ripening.Coarsening is associated with a power-law growth of typical domain size λ, which scales as λ ∼ t 1/3 with time t. 1 In the case of fluid-fluid phase separation, coarsening can also occur through advection with hydrodynamic flows that are driven by the interface tension between adjacent domains.Such hydrodynamic flows can speed up domain coarsening, which can lead to e.g. a linear scaling λ ∼ t. 1 However, the precise interplay between diffusive and hydrodynamic processes can be quite complex, leading to a rich phenomenology. 4ile phase separation has long been employed for industrial processes, it has also been identified as a possible key organizing mechanism in biological systems, from the sub-cellular 5 to the ecological 6 scale.For example, the interior of biological cells is organized in part through fluid compartments that can emerge from phase separation. 5,7,8Phase separation has moreover been used as effective model to describe the formation of multi-cellular aggregates. 91][12][13] This was motivated by experiments showing that in many animal species early differentiated cell populations unmix. 14,15ase separation of biological tissues may also help to understand an important question in the field of developmental biology: How does the vertebrate axis form during animal development?][18] Such gastruloids initially consist of spherical aggregates of undifferentiated stem cells (Fig. 1a).Even under isotropic boundary conditions, these aggregates then reliably 19 form polarized structures with different cell populations (Fig. 1b). 20This axis later defines the axis of the vertebral column. 17,21However, how exactly this axis initially emerges has remained unclear.
As a first step towards building a better physical understanding of such cellular aggregates, we study here the passive dynamics of a ternary phase separating system (Fig. 1c-e), where a droplet made out of two fluid phases (1 and 2) is suspended in a third fluid phase (3).Note that such multi-component emulsion droplets [22][23][24] are also important in many other contexts, including the food and pharmaceutical industries. 25There is a limited number of possibilities for the thermodynamic equilibrium state of such a ternary system, depending on the ratios of the interface tension σ 12 to the surface tensions σ 13 and σ 23 (Fig. 1f).Here, we focus on the case where the equilibrium state is a Janus droplet (Fig. 1d), i.e.where σ 12 < σ 13 + σ 23 and |σ 13 − σ 23 | < σ 12 (white region in Fig. 1f).While the equilibrium state is well known, we are interested in the dynamics that leads to this state, and how it depends on the interplay between diffusion and hydrodynamic flows.We start from an initial state where fluids 1 and 2 are almost homogeneously mixed inside of a spherical droplet (Fig. 1c).Thus, during the phase separation process, the system needs to spontaneously break rotational symmetry to reach the polarized equilibrium state (Fig. 1d).
Here, we numerically study the dynamics of this ternary system using 2D hybrid finite-volume Lattice-Boltzmann simulations of coupled Cahn-Hilliard and Stokes equations.We study the system's dynamics in terms of a Peclet number P e that characterizes the magnitude of hydrodynamic advection as compared to diffusive effects.We find that while increasing P e generally speeds up the phase separation process, intermediate P e give rise to longlived "croissant" states that slow down the relaxation to the final polarized state.These results are not fundamentally changed by including an asymmetry between the surface tensions σ 13 and σ 23 .Taken together, our work demonstrates the importance of hydrodynamic flows in deformable multi-phase droplets.

Free energy
We describe a 2D ternary system where the three components have area fractions C 1 , C 2 , and C 3 , which sum up to one everywhere: The phase-separating behavior of these three components is described by the following free energy: 26 ( The integral is over the total system area Ω.For each component j ∈ {1, 2, 3}, the free energy thus includes a quartic double-well potential with local minima at C j = 0 and C j = 1, and an energy density scale κ j .Moreover, gradient terms with coefficients κ j ensure finite widths of all interfaces.For any binary interface between phases of dominating components j and k, i.e.where the third area fraction vanishes, interface width α jk and interface tension σ jk take on the well-known values: 26 In the following, we will use the same width α = α jk for all interfaces jk ∈ {12, 13, 23}.

Dynamic equations
To define the dynamics of the area fractions C j such that Eq. ( 1) remains fulfilled, we follow Semprebon et al. 26 by introducing phase fields φ and ψ as Using these definitions together with Eq. ( 1), the free energy in Eq. ( 2) can be expressed in terms of φ and ψ instead of the C j .We then choose Cahn-Hilliard dynamics for both φ and ψ: where t is time, u u u is the local fluid velocity, Γ φ and Γ ψ are mobility parameters, and µ φ = δF/δφ and µ ψ = δF/δψ are the chemical potentials of φ and ψ, respectively.The factor [1 − ψ] in Eq. (5a) ensures that fluids 1 and 2 do not diffuse across fluid 3.This is motivated by the biological background -we exclude any dissociation and re-association of cells (components 1 and 2) from/to the aggregate.
To account for hydrodynamic effects, we describe all phases of the system as incompressible Newtonian fluids, using the incompressible Stokes' equations: Here, Π is the fluid pressure, η is the fluid viscosity, S S S is the shearrate tensor S S S = (∇u u u + (∇u u u) T )/2, and the last two terms correspond to the forces induced by the scalar fields. 1 We focus here on overdamped dynamics, as this is the relevant limit for ∼ 100µmscale cellular aggregates with shear rates on the ∼ day −1 scale.
Because the effective viscosity of embryonic cell aggregates is typically many orders of magnitude above that of water, 27 we choose a ψ-dependent viscosity η(ψ) = (1 − ψ)η i + ψηo, where η i ηo.

Initial and boundary conditions
We initialize the 2D system in a configuration schematized in figure 1c,e.The system initially consists of a circular droplet with diameter L, composed of a mixture of fluids 1 and 2 with ψ = 0, which is surrounded by fluid 3 with ψ = 1 and φ = 0.The inside of the droplet is initialized as φ = φ + ∆φ IC ξ, where φ and ∆φ IC are scalar parameters, and ξ is zero-mean, delta-correlated Gaussian white noise.The parameter φ corresponds to the conserved area average of φ, defined as φ = Ω (1 − ψ)φ dA/ Ω (1 − ψ) dA.We simulate the system in a squared periodic box with side length L t .Parameter values are listed in table 1.

Equilibrium state
The equilibrium state depends on four dimensionless parameters.These are first, the difference φ in area fractions between fluids 1 and 2.Here we focus exclusively on the case φ = 0, i.e. each fluid occupies half of the droplet area.Second, we set the ratio between droplet size L and interface width α to L/α = 64.
Finally, the equilibrium state also depends on the three tensions σ jk , for which we choose the following two dimensionless ratios: The parameter Ca is akin to a Capillary number -this ratio be- tween the inner and outer surface tensions controls the deformability of the droplet due to capillary effects.Unless stated otherwise, we set Ca = 1.The parameter ∆σ quantifies the surface tension difference between fluids 1 and 2 with respect to fluid 3.
In the following, ∆σ is varied over the range 0 ≤ ∆σ < 1.We are thus in the regime where the equilibrium state is a dipole structure (Fig. 1f).

Relaxation dynamics
The four phenomenological coefficients Γ φ , Γ ψ , η i , and ηo affect only the relaxation dynamics.Absorbing one to rescale time, we obtain three dimensionless parameters.First, we define a Peclet number that compares advective fluxes due to hydrodynamic flows to the diffusive fluxes of the φ-field: This number is defined on the scale of the droplet size L; it corresponds to the ratio P e = τ D /τ A between a diffusive time scale , and an advective time scale τ A = L/v f with the typical flow velocity scale v f = σ 12 /η i .Increasing the Peclet number corresponds to increasing the influence of advection.The limit case P e = 0 is simulated by ignoring the advection term in equations (5a)-(5b).In this case, we set Γ φ /Γ ψ = 3.In this article, we will vary P e to probe how the interplay of diffusion and hydrodynamic flows affects the phase separation dynamics.
Second, we define a Peclet number with respect to the ψ field, and with respect to the interface width α, as P e ψ = α 2 /(Γ ψ η i ).
Here, we fix this number to a relatively small value of P e ψ = 15/32 to ensure that the outer droplet interface remains stable during the simulations.
Third, we fix the viscosity ratio between droplet and the surrounding fluid to the value η i /ηo = 20.We have checked that changes of this ratio do not strongly affect the simulation results as long as η i /ηo 10.

Numerical implementation
We solve the dynamical equations (5a)-(6b) using a hybrid finitevolume Lattice Boltzmann approach.Similar methods have been used in earlier publications, [28][29][30] using finite-difference schemes for the phase field dynamics.Here, we use a finite-volume scheme instead to ensure exact conservation of the scalar fields φ and ψ. 31 We discretize space using a Cartesian square grid with spacing ∆x = α.For the advection terms in Eqs.central second-order scheme.We integrate over time using a firstorder explicit Euler scheme.We solve the incompressible Stokes equations (6a) and (6b) by including inertial terms with a small Reynolds number, which we integrate using a two-relaxation-time Lattice Boltzmann method. 32The time step for both finite-volume and Lattice Boltzmann parts is ∆t = (∆x/v f )/640 with the flow velocity scale v f = σ 12 /η i .Details on the numerical method are provided in the Supplementary Information.

Results
In the following, we examine how the droplet polarization dynamics depends (i) on the magnitude of advection as quantified by the Peclet number P e (sections 3.1 and 3.2) and (ii) on surface tension asymmetry, ∆σ = 0 (section 3.3).

Advection speeds up the polarization process
To study the role of advection for the phase separation process, we run simulations with symmetric surface tensions ∆σ = 0 and varying Peclet number P e (Fig. 2, Movies S1-S3).We find well-known features of spinodal decomposition, including an initial decomposition phase where |φ| generally grows until it saturates at ≈ 1, and a subsequent coarsening phase.However, Fig. 2 also suggests (1) that increasing the Peclet number speeds up the phase separation process (discussed below in this section),1 and (2) that intermediate Peclet numbers give rise to long-lived strongly elongated droplets (discussed in the next section, 3.2).
To quantify in how far increasing the Peclet number speeds up the droplet polarization process, we measure the time evolution of a dipole moment, which we define as Here, c c c is the barycenter of the droplet, defined as We study the normalized magnitude p = |P P P |/P 0 , where P 0 = L 3 /6 is a reference dipole moment, which corresponds to a circular droplet of diameter L split in two semicircular regions with φ = 1 and φ = −1.We find that indeed, the dipole moment p grows faster for higher Peclet number (Fig. 3a).We moreover define a polarization time tp of the system as the time when the polarization first reaches p ≥ 0.9peq, where peq is the dipole moment of the equilibrium state.Consistent with known results from the literature, 1 we find that advection speeds up the phase separation process (Fig. 3b).However, we note this effect only for P e 100, while tp remains approximately constant for Peclet numbers below 100.To better understand this P e dependence of the polarization process, we now discuss the different phases of the process separately.
Initial decomposition.We quantify the phase field amplitude |φ| during the initial decomposition phase, which lasts until the time t d ≈ 0.02τ D .To this end, we use the spatial root mean square φ (RMS) of φ, defined as φ2 = Ω (1 − ψ)φ 2 dA/ Ω (1 − ψ) dA.Fig. 3c shows that φ initially increases exponentially, and reaches a saturation at approximately t d .Except for very large P e, the behavior during this initial phase does not depend very strongly on the Peclet number.
Coarsening.We find that the subsequent coarsening phase is strongly affected by the Peclet number.This is apparent in Fig. 3d, where we plot the time evolution of the typical domain size λ (exact definition in SI).For P e = 0, the domain size λ coarsens according to an approximately constant power law with an exponent close to 1/3.In 2D, this is indeed the predicted exponent for systems that coarsen exclusively through the motion and deformation of domain boundaries, 1 even though we also observe occasional domain evaporation (Fig. 2).
For very large Peclet number, P e 10 3 , we observe a coarsening exponent close to one (Fig. 3d), which corresponds to the expected result for advection-dominated coarsening.For such large P e numbers, advection even induces noticeable coarsening already during the initial decomposition phase (Fig. 3d).This early coarsening is also the likely reason for the slowing down of the initial decomposition for very large P e (Fig. 3c), since diffusive unmixing needs to act over larger length scales in this regime.
The difference in the coarsening behavior between small and large Peclet number is also apparent when studying the behavior of the number of phase domains N D (Fig. 3e).We computationally define a phase domain as a connected region with constant sign of φ.The number of domains N D decreases during the coarsening process until it reaches a minimum of two, which is attained in the equilibrium state.Indeed, N D decreases more rapidly for larger Peclet number Fig. 3e.This confirms that advection speeds up the coarsening process, as N D and λ are strongly correlated (Fig. 3f).
For intermediate Peclet number, we initially observe a relatively fast coarsening (Fig. 3d).However, coarsening subsequently significantly slows down, which is also the reason why the polarization time tp decreases only for very large P e (compare Figs. 3b,d).What could be responsible for this slow down?

Elongated, striped droplets at intermediate P e
We wondered whether the elongated, striped droplets that emerge for intermediate P e (Fig. 2) could be related to the slow down of the coarsening at intermediate P e numbers.To test this hypothesis, we quantified the time evolution of droplet elongation E. We define E by expressing the droplet shape in polar coordinates and computing the second angular Fourier mode (Fig. 4a, details in SI).It is normalized such that E corresponds to the maximal radius variation of the second Fourier mode relative to the initial droplet radius r 0 = L/2.
We find that for P e = 0 and very large P e, the droplet elongation E monotonically increases with time until E ≈ 0.2 (Fig. 4b as high as E ≈ 0.3 before decreasing back to the final value of ≈ 0.2.Moreover, this regime of transient elongation is correlated in time to the phase of slow coarsening (compare Figs. 3d and  4b).This suggests that the elongated droplet structures could indeed be related to the observed coarsening slow down for intermediate P e.
Because the temporal coarsening dynamics differs when varying the P e number, we also plot the droplet elongation over the number of domains N D .This allows to compare droplet elongation E at similar "stages" of the coarsening process across different P e numbers.We find that at the same N D , elongation can be up to two orders of magnitude larger for intermediate P e (Fig. 4c).
We also observe in Fig. 2 that elongated droplets show striped domain patterns.Such a striped arrangement could indeed help to explain the observed coarsening slow down even for intermediate P e.A striped pattern of domains could be close to a mechanical equilibrium where only little hydrodynamic flows occur.In this case, coarsening would occur almost entirely through diffusion in an Ostwald-ripening-like process: Due to an increased Laplace pressure in small stripes as compared to larger stripes of the same color, the larger domain would slowly grow at the expense of the smaller domain through diffusive fluxes across the domains of opposite color.
To test these ideas and to better understand how elongated droplet shapes emerge, we quantify the nematic order of inner domain boundaries (i.e. between droplet domains 1 and 2) by an order parameter Q (Fig. 4d).To this end, we first introduce a tensor field A αβ that characterizes the local orientation of φ interfaces (Fig. 4d): Here, Greek letters are dimension indices and we use the nor- We then introduce the tensor Q αβ as the spatially averaged, symmetric, traceless part of A αβ : Here, Einstein notation is used and the brackets • denote spatial averaging over the system.We define the magnitude of the nematic order parameter as When the domain interfaces are randomly oriented, nematic order vanishes, Q = 0, while when all domain interfaces are straight and perfectly aligned, Q = 1.
Generally, the nematic interface alignment Q increases over time, and it does so more quickly for large P e numbers (Fig. 4e).Note that the final equilibrium state has Q ≈ 1, because it consists of a single straight internal interface.Notably, if we plot Q over N D , we find no strong difference across P e numbers (Fig. 4f).Indeed, Fig. 2 also show striped domain patterns not only for intermediate P e, but also for P e = 0.But if the domain pattern is also striped for P e = 0, then why do we find elongated droplets only for intermediate P e (Fig. 4g)?
A possible reason for why P e > 0 droplets are more elongated is that hydrodynamic flows would drive the system closer to mechanical equilibrium, where in mechanical equilibrium the droplets might be elongated.To test these ideas, we first analytically compute the mechanical equilibrium state for a striped droplet.For simplicity, we focus on the limiting case of an infinite chain of equally sized domains (Fig. 4h).In mechanical equilibrium with ∆σ = 0, interfaces between phases 1 and 2 are parallel straight lines, and surfaces to the external fluid 3 are circular arcs.One can furthermore show that the aspect ratio of a single domain in mechanical equilibrium is h/w = 2/Ca (derivation in SI).Thus, in a crude approximation, a droplet with N D domains in mechanical equilibrium would have aspect ratio To test whether droplets for intermediate P e approach the mechanical equilibrium state, we quantified the domain aspect ratio h/w in our simulations (details in SI).We plot h/w for Ca = 1.5 depending on P e in Fig. 4i and compare to the theoretical solution of the infinite chain (dashed line).Our results suggest that the droplets do indeed approach mechanical equilibrium as P e increases.
To further test these ideas, we compared domain aspect ratios also across different values of Ca (Fig. 4i inset).We found that for smaller Ca, the measured aspect ratio started to become smaller than the theoretical prediction.This makes sense since for smaller Ca, the overall droplet aspect ratio also becomes smaller and so the simulated droplets deviate more strongly from the theoretical picture of a chain of equally sized domains.Indeed, for Ca → 0, the infinite-chain prediction would be h/w = 2/Ca → ∞.However, in this limit, there is no inner interface tension, σ 12 = 0, and the droplet shape is spherical.When dividing a spherical droplet into N D domains of equal width, the domain aspect ratio is h/w ≤ N D , which corresponds to an approximate upper bound for h/w.
Taken together, our results show that for intermediate P e, transient droplet states form, which are striped and elongated.They are close to mechanical equilibrium and coarsen only through diffusive fluxes, which makes them long-lived.

Effect of a difference between the surface tensions of the two droplet phases
We next study the effect of different surface tensions between the two droplet phases, ∆σ = 0. To this end, we carried out simulations with varying 0 ≤ ∆σ < 1 for the Peclet numbers 0, 100, and 10 4 (Fig. 5, Movies S1-S9).
We first examined how varying the surface tension difference ∆σ would change the time tp required to reach the polarized equilibrium state (Fig. 6a).For P e = 0, only simulations with ∆σ = 0 and ∆σ = 0.4 reached the polar equilibrium state within the maximal simulation time of tmax = 50.For P e = 100 and P e = 10 4 , the equilibrium state has always been reached after a finite time tp, but no clear trend in the dependency on ∆σ is apparent.To better understand these results, we again separately discuss several phases of the process.

Initial decomposition.
For ∆σ > 0, the initial decomposition into the two droplet phases proceeds from the droplet boundary inwards in concentric circles (Fig. 6b).This phenomenon is called composition wave in the literature: 33,34 For ∆σ > 0, fluid 2 (orange) has a lower surface tension, and so a stripe with an abundance of fluid 2 first appears close to the surface.This supplants fluid 1 (blue) in this region, pushing it further inwards.The slight abundance of fluid 1 next to the outer fluid-2 stripe attracts more of fluid 1, supplanting more of fluid 2, some of which is pushed further inwards, and so on.This initial decomposition proceeds until the whole droplet displays a rotationally symmetric striped pattern, where the stripes reach saturation |φ| ≈ 1 again around Coarsening.The rotationally symmetric pattern coarsens through progressive broadening of the stripes and occasional annihilation of the innermost phase domain (Figs. 5, 6c).For most simulations with P e = 0, the system most often remained in a state with two concentric phase domains until the maximal simulation time (e.g.Fig. 5, for P e = 0, ∆σ = 0.8).However, for finite Peclet number, this rotational symmetry got broken; for P e = 100 during the coarsening phase, and for P e = 10 4 already during the initial decomposition phase.time when the nematic order of the inner domain boundaries, Q, surpassed a value of 0.01.For P e = 100, we found that t b increased monotonically with ∆σ (Fig. 6d & inset).
Transient flows.The symmetry breaking usually triggered complex transient hydrodynamic flows that were qualitatively similar in different parts of the parameter space, and reminiscent of Marangoni flows within a droplet (Fig. 6f). 35,36For P e = 100, these flows typically resulted in a croissant state.
Slowly relaxing "croissant" state.We wondered whether also the time te it takes for the croissant state to emerge shows a clear dependence on ∆σ.We defined te as the first time when the droplet elongation E reached 90% of its maximum value during any given simulation.This captures both the approximate time of maximum elongation, in cases where a croissant state was created, or the time when E first reached the final plateau value, in cases where the transient flows directly created a 2-phase Janus droplet.We find that te does indeed increase monotonically with ∆σ (Fig. 6e & inset).Thus, while the time tp for the whole process varies nonmonotonically with ∆σ, the time te when the croissant state appears increases monotonically, and is generally much smaller than tp (compare Figs. 6a and e).This means that the variation in tp is due to a variation in the life time of the croissant states.
In the previous section, we concluded that for intermediate P e, these states coarsen only due to an Ostwald-ripening-like process, where smaller domains shrink at the expense of larger ones.This suggests that the variation in tp may just be due to the way the two droplet phases happen to be distributed into different domains by the transient flows that follow the symmetry breaking.Indeed, we find for instance that for P e = 100, the simulations with ∆σ = 0.4 need longest to relax to the 2-domain equilibrium state (Fig. 6a), and these are also the simulations where two approximately equally-sized orange domains form (Fig. 5), leading to an initially slow diffusion between the two domains.
Finally, in simulations with P e = 10 4 , a first symmetry breaking event occurred at a time around 0.01.The ensuing transient flows often resulted in a second approximately rotationally symmetric state with 2 domains at a time around 0.02.This state only later resolved into the polar equilibrium state.We believe this second, 2-domain rotationally symmetric state forms in this case because the first breaking of rotational symmetry occurs before the initial decomposition is finished, i.e. when |φ| is significantly smaller than one within the domains.As a consequence, the interface tension between both domains, which scales as ∼ φ 2 , is smaller, and thus the effective value of ∆σ = (σ 13 − σ 23 )/σ 12 is larger than its nominal value.Hence, we believe that we find the rotationally symmetric 2-domain state as transient state, because the expected equilibrium state is the rotationally symmetric one as long as the effective ∆σ is larger than one (Fig. 1f).

Discussion
In this article, we discussed phase separation in a two-phase deformable droplet.While the thermodynamic equilibrium state for such a system is well known, we studied the route to reach equilibrium.In particular, we focused on the interplay between spinodal decomposition and advection with hydrodynamic flows created by interface tensions.While this interplay has been studied before, 4 it was completely unknown how it plays out in a finite, deformable system.We found that advection can speed up the coarsening process as reported earlier. 1 We moreover found that for intermediate P e numbers, elongated, striped droplet states emerge, which we call "croissant" states.The stripes correspond to a nematic alignment of inner domain boundaries, which emerges for low and intermediate P e.For intermediate P e, the striped droplets elongate as they approach mechanical equilibrium.These states are long-lived, because they coarsen almost exclusively via diffusion, through an Ostwald-ripening-like process.We did not observe such "croissant" states for very large P e numbers, where coarsening is entirely dominated by advection, which leads to coalescence of the phase domains before a striped pattern can develop (Movie S3).
We also studied the effect of an asymmetry between the surface tensions of the two droplet phases, and we found that it changes in particular the beginning of the phase separation process.The emergence of composition waves initially creates a rotationally symmetric system, which then starts to coarsen diffusively.Breaking of this rotationally symmetric state triggers transient flows that are reminiscent of Marangoni flows in droplets. 35,36These transient flows typically gave rise to a long-lived croissant state.
Taken together, our results demonstrate that advection can play an important role in the coarsening dynamics of deformable twophase droplets.

Symmetry breaking in stem cell aggregates
][39] Aggregates of stem cells that mimic early vertebrate axis formation are known to develop a polar structure (Fig. 1a,b), 16 reminiscent of two-phase Janus droplets in thermodynamic equilibrium (Fig. 1d), where the different phases correspond to different cell types.While current studies of such stem cell colonies focus almost exclusively on the biochemical dynamics, 40,41 our work here demonstrates that tissue flow and advection could play an important role as well.Indeed, recent experiments have identified significant large-scale tissue flows in stem cell aggregates. 20s a first rough estimate, the Peclet number in mouse embryonic stem cell aggregates is on the order of P e ∼ Lσ/Dη ∼ 200.This is based on an aggregate size of L ∼ 200 µm, 20 a tensionto-viscosity ratio of σ/η ∼ 5 µm/h, 42 and a diffusion constant of order D ∼ 5 µm 2 /h (unpublished preliminary data, lab of Pierre-François Lenne).In this P e regime, we observe in our simulations the emergence of the "croissant" states.Indeed, these states appear visually reminiscent of somites, which have also been observed to form in gastruloids, 21 while an interface tension has been reported at somite-somite boundaries in vivo. 43However, first, the mechanism of somite formation in vertebrates is generally believed to crucially involve the complex interplay of biochemical signals in a clock-and-wavefront model 44 (even though purely mechanical models are discussed as well 45 ).Second, different from our model, gastruloids exhibit a polar organization before somites form. 21This raises the question why before polarization there are usually no "croissant" states observed in experiments.
Stem cell cultures are of course much more complex than the model we studied here, and future refined studies should take additional effects into account.First, experiments indicate that the dynamics of some biological signaling molecules and transcription factors is non-conserving in these systems. 20,40,41To take this into account in our model, a source term needs to be added in the φ equation Eq. (5a).Similarly, we studied here a droplet of constant volume, while the stem cell colonies are typically growing.To implement this, source terms need to be added to the ψ dynamics Eq. (5b) and in the incompressibility condition Eq. (6a).
Second, we focused here on a 2D model, while the stem cell systems motivating this study are 3D aggregates. 16,20While there are likely several commonalities between the 2D and 3D cases, we expect differences for instance due to Laplace pressure effects. 34,46,47Such effects arise only in 3D, where in tube-shaped phase domains, interface tension creates a Laplace pressure that depends on the tube radius.This pressure can lead to necking and breakup of connected domains in a Plateau-Rayleigh instability, 46 but it can also lead to a pumping effect that makes the surface layer of the droplet grow more rapidly. 34It will be interesting to dissect the consequences of these effects for a deformable 3D droplet in future work.
Third, the cell aggregates are active systems, and so activity may need to be included for an effective description of symmetry breaking in the aggregates. 48,49ourth, additional hydrodynamic fields may be required to fully account for some features of the stem cell aggregate dynamics.In particular, the inclusion of more scalar fields in addition to φ may be required, 40,41 or of a polar field as experiments on fish aggregates have pointed to a role of polarity proteins. 50,51inally, the real system dynamics is likely influenced by several types of noise.Moreover, a φ-dependence of mechanical properties such as viscosity or local volume growth may be required to explain some of the experimental observations. 52

Fig. 1
Fig. 1 Phase separation and symmetry breaking in deformable multicomponent droplets.(a,b) Microscopy images of cell aggregates mimicking early vertebrate development, taken 1 day apart from each other.Brighter regions indicate the early mesoderm marker T/Brachyury (images taken from Ref. 20 ).(c-e) Ternary fluid mixture model: a droplet, immersed in fluid 3 and initially composed of a mixture of fluids 1 and 2 (c), undergoes phase-separation until reaching the thermodynamic equilibrium state (d) characterized by a surface tension balance between all interfaces.(e) Schematic representation of the physical configuration.(f) Dependence of the thermodynamic equilibrium state on a Capillary number Ca = 2σ 12 /(σ 13 + σ 23 ) and the relative surface tension difference ∆σ = (σ 13 − σ 23 )/σ 12 .The hatched region is unphysical, because at least one of the tensions σ ij would need to be negative there.

1 - 11 J
o u r n a l N a me , [ y e a r ] , [ v o l .] ,

Fig. 2
Fig.2Advection affects the coarsening dynamics and the time to reach the polarized equilibrium state.For each value of the Peclet number, a time series of snapshots is shown, representing the φ field from −1 (blue) to 1 (orange).Here we set ∆σ = 0, i.e. the outer surface tension is independent of φ.

Fig. 3
Fig. 3 Advection speeds up droplet polarization by accelerating the coarsening phase.(a) Dipole moment p depending on time for different P e.(b) Polarization time tp as a function of the Peclet number; error bars indicate the standard error of the mean.(c) Spatial root-meansquare of φ depending on time.The black line shows the initial growth with rate ω predicted by linear stability analysis, φ ∼ e ωt .(d) Domain size λ as defined in the SI depending on time.(e) Number of phase domains N D depending on time.(f) Domain size λ as a function of the number of phase domains N D .Legends in panels c-f same as in panel a.In panels c-e, the vertical dashed line indicates the approximate end time t d ≈ 0.02τ D of the initial decomposition phase.The quantities on all ordinate axes are averaged over 25 simulation runs.

4 10 0 10 1 10 2 10 3 10 4 N D = 5 NFig. 4
Fig. 4 Transient "croissant" states of the droplets, which show an elongated and striped organization of the phase domains.(a) Illustrative example of droplet surface deformation, at P e = 100 and t = 0.1, and resulting angular Fourier spectrum used to define the droplet elongation E. (b) Evolution of the droplet elongation as a function of time.(c) Transient droplet elongation as a function of the Peclet number at fixed numbers of phase domains, N D .(d) Illustrative nematic field based on the φ-field gradient, at P e = 100 and t = 0.1.(e,f) Evolution of the global nematic order as a function of (e) time and (f) N D , for different Peclet numbers.(g) Comparison for two snapshots for P e = 0 and P e = 10 with N D = 5.(h) Limiting case of infinite chain of equally-sized phase domains in mechanical equilibrium, where h/w = 2/Ca (see SI). (i) Phase domain aspect ratio h/w as of function of P e for Ca = 1.5.(i inset) Aspect ratio over Ca for P e = 10.The aspect ratio is averaged for time points with N D = 5 over all domains having exactly two neighboring domains (i.e.all non-branching inner domains).The quantities in all panels are averaged over a series of 25 simulations.

7 Δσ=0Pe=10000Fig. 5
Fig. 5 Surface tension difference affects transient droplet morphologies and polarization time.For different values of ∆σ and P e, a time series of snapshots is shown, representing the φ field from −1 (blue) to 1 (orange).

Fig. 6 (
Fig. 6 (a) Polarization time as a function of ∆σ.(b) Composition wave during the initial decomposition phase, for P e = 100 and ∆σ = 0.4.(c) Domain size λ as a function of time.(d) Symmetry breaking time t b over ∆σ.The symmetry breaking time t b is the first time when the nematic order Q of inner domain boundaries increases above a value of 0.01.(d inset) Nematic order Q of the inner domain boundaries over time for different ∆σ.Legend same as in panel c.(e) Time te when a croissant state (or a Janus droplet) emerges over ∆σ.The time te is defined as the first time when droplet elongation E reaches 90% of its maximum value.(e inset) Droplet elongation E over time for different ∆σ.Legend same as in panel c.(f) Directional flow during symmetry breaking.In panels b-f: P e = 100.In each of the panels c-e incl.insets, the respective quantity on the ordinate axis is averaged over 25 simulation runs.

Table 1
Values of dimensionless parameters used in our simulations