Discontinuous switching of position of two coexisting phases

Here we investigate how the positions of a condensed phase can be controlled by using concentration gradients of a regulator that influences phase separation. We consider a mean field model of a ternary mixture where a concentration gradient of a regulator is imposed by an external potential. We show that novel first order phase transition exists at which the position of the condensed phase switches in a discontinuous manner. This mechanism could have implications for the spatial organization of biological cells and provides a control mechanism for droplets in microfluidic systems.

Phase separation of a mixture refers to the formation of a condensed phase that coexists with a dilute phase of lower concentration [1,2]. Such demixing is the result of a first order thermodynamic phase transition where the concentration difference between the phases changes discontinuously. It can be observed in many forms in everyday life, for example when oil is added to water. The occurrence of a transition from the homogeneous mixture to a system with coexisting phases can be controlled by temperature or by changing the composition of the mixture. Condensed phases are influenced by surfaces possibly causing wetting transitions [3][4][5]. Furthermore, phase separation can be affected by external forces such as gravity causing sedimentation. Here we study the equilibrium physics of the positioning of two condensed phases in inhomogeneous systems. We present a simplified model that provides the basic mechanism for the positioning at thermal equlibrium which can be further extended to non-equilibrium processes such as the kinetics of droplet formation and ripening. In our model phase separation of two components is subject to a concentration gradient of a regulator component where the gradient is generated by an external field. The regulator component affects demixing of the two components but does not phase separate itself. The system then relaxes to a spatially inhomogeneous thermodynamic equilibrium state with two coexisting phases positioned by the regulator gradient. The spatial distributions of the three concentration profiles at thermal equilibrium are determined by minimizing a mean field free energy functional.
We find that as a function of an interaction parameter the position of the condensed phase switches discontinuously from a position in the region of large regulator concentration (correlated state) to the region of low regulator concentration (anti-correlated). This switching R (x) (a) (b) of position corresponds to a novel, equilibrium first order phase transition at which an order parameter undergoes a jump ( Fig. 1(a,b)). In our equilibrium model for spatial regulation of phase separation we consider three components [10]: two components which can demix from each other, A and B, and a regulator R that interacts with these components. The regulator affects phase separation but does not demix from A and B. Demixing and interactions with the regulator are described by the Flory-Huggins free energy density for three components ( [11,12] and Supplemental Material [13], II): We consider the incompressible system in which the molecular volumes are equal to ν and The logarithmic contributions correspond to the mixing entropy, while the second line in Eq. (1) describes the molecular interactions between the components; χ ij is the interaction parameter between component i and j. The gradient terms represent contributions to the free energy associated with spatial inhomogeneities. They introduce two length scales, √ κ A and √ κ R . The regulator R is subject to an external field described by a position-dependent potential U (x). For simplicity we consider in the following a onedimensional system and choose a potential of the form where s > 0 characterizes the slope of the potential and its inverse is a third length scale in our model. Note that in the absence of A and for φ R φ B , φ R (x) attains a concentration profile that is linear in space with a slope s. We consider a finite system of size L and two type of boundary conditions: (i) Neumann boundary conditions, for all fields, where the primes denote spatial derivatives, and (ii) periodic boundaries with . The conditions (i) imply that there is no explicit energetic bias to wet or dewet the boundary, but the presence of the boundary enforces the slopes of the concentration profiles close to the boundary. In contrast, the periodic conditions (ii) allow to study the system in the absence of boundaries.
To calculate the equilibrium profiles φ A (x) and φ R (x), we minimize the free energy Due to particle number conservation, two constraints are imposed for the minimization: , whereφ i are the average volume fractions andφ B = 1 −φ A −φ R . Variation of the free energy Eq. (2) with the constraints of particle number conservation implies (i = A, R): where λ R and λ A are Lagrange multipliers, and the prime denotes a derivative with respect to x . The boundary terms vanish for both, Neumann and periodic boundary conditions.
Using the explicit form of the free energy density (Eq. (1)), the Euler-Lagrange equations can be derived (see Supplemental Material [13], I)). We solve these equations using a finite difference solver (bvp4c in MATLAB [14]). As control parameters we consider the three interaction parameters χ AR , χ AB and χ BR , the slope of the external potential s and the mean volume fraction of A-material,φ A . The mean regulator material is fixed toφ R = 0.02 in all presented studies. Moreover, we focus on the limit of strong phase segregation where the interfacial width is small compared to the system size, i.e. √ κ i L. In this limit, we verified that our results depend only weakly on the specific values of κ i . interaction parameter χ BR . F l and F r are the free energies of the correlated and anti-correlated stationary solution with respect to the regulator gradient, respectively ( Fig. 1(a,b)). Lines are dashed when solutions are metastable. At χ * BR , F l and F r intersect causing a kink corresponding to the solution of lowest free energy. This shows that the transition between correlation and anticorrelation is a discontinuous phase transition. (b) The order parameter ρ BR (Eq. (4)) jumps at χ * BR by a value of ∆ρ * BR . The transition point χ * BR does not depend on the slope of the regulator s, while χ * BR increases linearly with s (see Supplemental Material [13], VI) Parameters: χ AB = 4, χ AR = 1,φ A = 0.5,φ R = 0.02, κ R /L 2 = 7.63 · 10 −5 , κ A /L 2 = 6.10 · 10 −5 , κ/L 2 = 6.10 · 10 −5 , Ls = 0.99. For plotting, ν = L/256 was chosen. Solving the Euler-Lagrange equations with Neumann boundary conditions (i), we find two spatially inhomogeneous solutions for component A, which we denote φ l A (x) and φ r A (x), and the two corresponding solutions for the regulator component R, are denoted φ l R (x) and φ r R (x) (the profile of B follows from volume conservation). The phase separating material A is either accumulated close to the right boundary of the system (φ r R (x) and φ r A (x)) and correlated with the concentration of the regulator material ( Fig. 1(a)), or it is accumulated at the left (φ l R (x) and φ l A (x)) and anti-correlated with the regulator ( Fig. 1(b)). Upon varying the interaction parameters χ BR in Fig. 7(a,b), the free energies of the correlated and the anti-correlated states, Fig. 7(a)). At this point the lowest free energy exhibits a kink, which means that the system undergoes a discontinuous phase transition when switching from the spatially anti-correlated ('left') to the spatially correlated ('right') solution with respect to the regulator. A set of order parameters suitable to study this phase transition is where the squared normalization N 2 , where Θ(·) denotes the Heaviside step function. This normalization ensures that −1 < ρ ij < 1 and . The derivative of the free energy with respect to the interaction parameter χ ij generates the covariance between the spatially dependent fields φ i (x) and φ j (x). If the fields are spatially correlated, ρ ij > 0, and if they are anti-correlated, ρ ij < 0. For homogeneous fields with Varying the interaction parameter χ BR (Fig. 7(b)), the order parameters ρ BR and ρ AR jump at the threshold value χ * BR , while in the absence of a regulator gradient (s = 0), they change smoothly. The jump of both order parameters in the presence of a regulator gradient indicates that the spatial correlation of A and B to R changes abruptly, which is expected in case of a first order phase transition.
By means of the order parameter ρ BR (Eq. (4)) we can now discuss the phase diagrams as a function of the interaction parameters for different volume fractions of the demixing material,φ A . We find three regions ( Fig. 3(a-c)): A mixed region (M), where volume fraction profiles are only weakly inhomogeneous and no phase separation occur. In addition, there are two regions, (C) and (AC), where components A and B phase separate and A is spatially correlated or anti-correlated with the regulator R, respectively. There exists a triple point where all three states have the same free-energy. Forφ A = 1/2, the shape of the transition line between correlated and anti-correlated states is straight and χ * BR is independent of χ AB (Fig. 3(b)). Ifφ A is decreased, the region of the correlated state in the phase diagram grows.
In this case, the correlated state is favored, while for increasingφ A , the anti-correlated state is preferred. The transition line to the mixed states is horizontal forφ A = 1/2 ( Fig. 3(b)).
For both, larger and smallerφ A -values, it becomes curved and moves towards larger χ AB interaction parameters. This behavior can be qualitatively understood by the upshift of the demixing threshold χ AB onceφ A deviates from 1/2, as known for binary systems. Since the concentration of R is small here, this analogy provides a good approximation (φ R → 0 in Eq. (1)). Both trends explain the parabolic shape of the positions of the triple point in the phase diagrams whenφ A is varied (Fig. 3(d)).
The transition line in the phase diagrams between the correlated and anti-correlated solution as a function of the interaction parameters can be estimated analytically. In the absence of a regulator gradient (s = 0), the free energies of both solutions are the same for all interaction parameters for which phase separation occurs. In the presence of a regulator gradient, however, the free energies corresponding to the correlated and the anti-correlated solutions are unequal for most points in the phase diagram. The reason is that the external potential U (x) forces the regulator to form a gradient, and thus the interactions with the regulator lead to different free energies of the correlated and anti-correlated states. Only along the transition line between both states the free energies equal: This condition can be used to estimate the transition line for varying interaction parameters and the slope of the potential, s. To estimate ∆F we parametrize the profiles of the stationary solutions φ r,l A (x) and φ r,l R (x) using physical assumptions that are in agreement with our numerical results. First we idealize the already narrow interface of the demixed component φ A as sharp. Since the regulator is maintained by the external potential, we find φ r R (x) φ l R (x) close to the transition line. Thus we use the one profile, denoted as φ R (x), for both regulator states. In addition, we approximate the regulator profile as linear function with slope m, neglecting spatial non-linearities that can be seen in Fig. 1(a,b). The low volume fractions outside the condensed phase of the demixed binary A-B system are approximated as constant valuesφ out . The larger volume fraction (inside) shows a weakly linear profile ( Fig. 1(a,b)). For most parameters, the volume fraction inside the condensed phase can be well described as φ in (x) =φ in − φ R (x), whereφ in is the constant volume fraction inside the condensed phase of the binary A-B mixture (see Supplemental Material [13], V). The approximated profiles are: The conservation of A determines the domain sizes l,r of the phase separated region (see Supplemental Material [13], IV). To calculate ∆F (Eq. (5)), the free energy density (Eq. (1)) is integrated in the domain [0, L]. Using the approximated profiles (Eqs. (6)) we find where the value G depends only on the parameters of the simplified solutions (see Supplemental Material [13], IV). Consistently, ∆F = 0, if there is no regulator gradient (m = 0).
In presence of a regulator gradient, ∆F = 0 if χ * BR = χ AR , which defines the transition line between the correlated and anti-correlated solution obtained from the parametrized solutions Eqs. (6). This prediction is in very good agreement with our numerical results for φ A 1/2; see black lines in Fig. 3(a) and Fig. 4(a). By means of the ansatz given in Eqs. (6) we can also estimate how the jump of the order parameter ∆ρ * BR (definition see Fig. 7(b)) at the transition point depends on the model parameters. In particular we find that the estimated ∆ρ * BR as a function of the slope of the regulator (not shown) and the interaction parameter χ AB (Fig. 4(b)) almost perfectly describe the data obtained from the numerical minimization of the free energy. This shows that the proposed parametrization of the stationary solutions represents a consistent approximation. We conclude that the positioned and phase separated profiles possess a sharp interface and the volume fraction inside has a weak linear slope that is mainly determined by volume exclusion with the regulator.
The phase diagrams (Fig. 3) depend on the boundary conditions rising the question whether the boundary play a key role for the existence of the phase transition. To this end we considered a periodic system without boundaries. We find that the reported first order transition also exists for in the absence of boundaries (see Supplemental Material [13], III). Thus the transition is not induced by boundaries as for example in the case of wetting transitions [3][4][5].
The discontinuous switching of phase separation could be tested experimentally. A soluble salt of high magnetic susceptibility could be used to create and maintain concentration gradients via the application of an inhomogeneous magnetic field [15]. Phase separation in a regulator gradient could be observed by introducing components that phase separate in a salt dependent manner. In particular, a pre-formed droplet could be added to an existing regulator gradient or the regulator gradient is created after Ostwald-ripening is completed [16][17][18][19]. The phase transition could be triggered by changing the concentrations of the phase separating material, by changing the temperature or by adding additional components that influence the interaction parameters. The systems considered here could also be relevant for applications. As the composition of a condensed phase creates a distinct chemical environment, our work may provide a novel mechanism to control and switch chemical environments in microfluidic devices.
We would like to thank Martin Elstner and Omar Adame for fruitful and stimulating discussions. This project was supported by the Center for Advancing Electronics Dresden (cfAED). Christoph A. Weber thanks the German Research Foundation (DFG) for financial support. Samuel Krüger and Christoph A. Weber contributed equally to this work.

Euler-Lagrange Equations
From the variation of the free energy and the explicit form of the free energy density (main text, Eq.(1)), we find the following set of the Euler-Lagrange equations: Here, we defined χ = χ AR − χ AB − χ BR and rescaled length x → x L.
Penalty of spatial inhomogeneities in the ternary Flory-Huggins free energy density

Derivation using a mean field approximation
To show this relation, we start from the local mean-field free energy on the lattice and calculate the continuum limit of this free energy as shown in Ref. [20] for a binary system.
The local free energy density of the three component system is derived in [21,22] using a mean-field approximation: where ν is the molecular volume. The greek indices α and β indicate the positions on the lattice. The first line describes the entropy of the mixture. Each contribution is local. The second and third line contains the energetic part of the free energy. It describes the non-local interactions between neighboring lattice sites.
In the next steps we will perform the continuum limit. In case of the entropic contribution, we can simply replace φ i (α) → φ i (x). In case of the energetic contributions, we rearrange the terms leading to: Each contribution can be rewritten as We can identify the Flory Huggins interaction parameter as χ ij = 1 2 β J ij (α, β). In the continuum limit we can introduce the gradient of the volume fractions as (φ i (α) − φ i (β)) → a∇φ i . We finally obtain the free energy F = dx f with the free energy density given as where The parameters characterizing the penalty corresponding to spatial inhomogeneities are

Phenomenological derivation
In the Ginzburg-Landau free energy the penalties corresponding to spatial inhomogeneities are phenomenologically introduced based on symmetry considerations: whereκ i > 0 since spatial inhomogeneities are unfavored. Moreover, f 0 is the free energy density that only depends on the volume fractions φ i , i ∈ A, B, R. However, only two volume fraction fields are independent due to particle conservation and incompressibility, Here, κ A =κ A +κ B , κ R =κ R +κ B and κ =κ B .

Choice of the parameters κ i
In the presented studies, we have chosen κ A = κ for simplicity. Please note that the derivation presented in Sect. is based on a mean field approximation and therefore it should only serve as an estimate for the values κ i . We chose the values for the parameters κ A and κ R consistent with these estimates (see figure captions in the main text).

Discontinuous phase transition in a periodic domain
Here we discuss the results of the minimization of the free-energy (Eq. (3), main text) using periodic boundaries with φ i (0) = φ i (L) and φ i (0) = φ i (L). We find the same main results as for Neumann boundary conditions, namely the existence of a discontinuous phase transition. In the periodic domain, we also use a periodic external potential: The parameter ω is a phase shift. The value of the phase is chosen such that that the region of segregated A-material is placed at x = 0. The logarithmic form of the potential is chosen ensures that a sinus distribution of the regulator is obtained in the dilute limit. We find two stationary solutions of different spatial correlations with respect to the regulator. They switch at χ * BR by a discontinuous phase transition ( Fig. 5(a-c)). Therefore, a boundary of the system is not a necessary requirement for the emergence of the discontinuous phase transition discussed in our manuscript.

Estimate of ∆F
The free energy difference between the two stationary solutions, ∆F , results from integration over the domain [0, L] using the simplified solutions (Eqs. (7), main text): where and Here we substituted the interaction parameter between B and A by χ * AB , and truncated at O(∆ ) with ∆ = r − l . Consistently, ∆F = 0, if there is no regulator gradient (m = 0), and when phase separation is absent ( l = r = 0, L). G depends only of the parameters of the simplified solutions (see Eqs. (7), main text). The jump of the order parameter at the transition point linearly increases with the slope of the gradient s. The slope of this linear dependence is influenced byφ A . Fixed parameters: χ AB = 4, χ AR = 1,φ R = 0.02, κ R /L 2 = 7.63 · 10 −5 , κ A /L 2 = 6.10 · 10 −5 , κ/L 2 = 6.10 · 10 −5 , ν = L/256.

Regulator Peak at the Interface
The numerically obtained regulator profiles show a significant peak at the interface between the A-rich and the B-rich phase (see main text, Fig. 1(a,b)). The emergence of the regulator peak can be understood by entropic and energetic considerations of the free energy.
For large and positive χ AR and χ BR (corresponding to a repulsive tendency with respect to the regulator), the energy of the system decreases as regulator accumulates at the interface.
Moreover, the entropy decreases as the composition of the interfacial region of all three components is closer to a well-mixed state. The amount of regulator material that is accumulated at the interface is strongly influenced by the κ i -parameters; see Fig. 8. In Fig. 9 left, the peak area is shown for varying κ i -parameters. For simplicity, we chose κ A = κ R = κ. The peak area vanishes as the κ iparameters approach zero. This behavior is expected since these parameters set the size of the interface between the phase separated phases. In this limit, the estimates for the phase boundaries based on the approximate solution (main text, Eq. (6)) are valid.
However, the peak height and thereby the existence of the peak is approximately independent of κ i (Fig. 9 right). This indicates that the existence of the peak may depend on the interaction parameters for example. Since we also observed that the peak is more pronounced at the transition line between anti-correlated state and correlated state, we investigated the energetic influence on the peak height along the transition line. As derived in the main text, the transition line is governed by the condition χ AR = χ BR forφ A = 0.5.