Phase separation and morphology formation in interacting ternary mixtures under evaporation: Well-posedness and numerical simulation of a non-local evolution system

We study a nonlinear coupled parabolic system with non-local drift terms modeling at the continuum level the inter-species interaction within a ternary mixture that allows the evaporation of one of the species. In the absence of evaporation, the proposed system coincides with the hydrodynamic limit of a stochastic interacting particle system of Blume-Capel-type driven by the Kawasaki dynamics. Similar governing dynamics are found in models used to study morphology formation in the design of organic solar cells, thin adhesive bands, and other applications. We investigate the well-posedness of the target system and present preliminary numerical simulations which incorporate evaporation into the model. We employ a finite volumes scheme to construct approximations of the weak solution and illustrate how the evaporation process can affect the shape and connectivity of the evolving-in-time morphologies.


Introduction
In this paper, we study the solvability and the numerical simulation of a nonlinear coupled parabolic system with a specific shape of the non-linear non-local drift, which is derived in [24] as hydrodynamic limit of the Blume-Capel model endowed with Kawasaki dynamics describing the interaction of ternary mixture of particles.Such mixtures typically consist of two different solutes mixed within a solvent and allow for phase separation.The situation is rather typical in modern materials science when thin films are involved; see e.g., Section 1.2 for a discussion of the application of such scenario to producing morphologies as needed for organic solar cells.In the applications discussed below, the solvent particle evaporates and plays a crucial role in the morphologies produced by the separated phases.The presence of such a non-equilibrium process makes this scenario different compared to the more classical phase separation settings arising in crystal growth, metallurgy, and so on, normally treated by Allen-Cahntype or Cahn-Hilliard-type equations.Specifically, by controlling the evaporation mechanism a local equilibrium is frozen, i.e., one has the possibility to actively select specific morphologies as the end configuration.We refer the reader to the recent works [6,28,34], where this setting has been explored using Monte Carlo simulations for suitable lattice-based models.Although we are driven by these two concrete applications, our work can be seen as a continuation of the efforts started earlier by Lebowitz and Giacomin [14,15].
For a good grip on the design of morphologies out of ad hoc polymer-polymer-solvent interactions, it is crucial to have sufficiently rich models posed at the continuum level that inherit, from the interacting particle systems, the ability to produce physically meaningful phase separation.This would set the stage for process and shape optimization operations to reach optimal macroscopic transport and reaction properties, mediated by optimal morphology shapes.However, in spite of huge efforts in the statistical mechanics and applied mathematics communities, rigorous derivations of the corresponding hydrodynamic limit equations for mixtures of interacting particle systems with evaporation are currently out of reach.The main technical difficulties in performing rigorously the hydrodynamic limit process seem to arise because of the evaporation component of the process.To allow for the evaporation of the solvent, we modify the limit system with a linear evaporation term to mimic the release of the solvent.The problem setting is relevant, with a rich phenomenology, for the presence of three coexisting phases undergoing the separation process.

A model without evaporation
Within this framework, the main interest is on the solvability of the following system of two coupled non-local nonlinear parabolic equations Here, t and x represent the time and space variable, respectively, m represents the average spin density (also called magnetization), and ϕ represents the solute volume concentration (we expand upon the meaning of these terms in a later section).Additionally, Ω ⊂ R d is a cube with spatially periodic boundary conditions (i.e., homeomorphic to the unit torus), T is some positive time, β > 0 is a constant representing the inverse temperature, J ∈ C 2 + ( Ω) is a symmetric compactly supported potential such that Ω J(r)dr = 1.We note that in the analysis to follow, such regularity on J is unnecessary.Indeed, the proofs below readily hold for J ∈ W 1,∞ (R d ) and small adjustments can be made for to include J ∈ W 1,2 (R d ) and possibly J ∈ W 1,1 (R d ); however, we maintain the C 2 regularity to stay true to the Kac interpretation of the potential.For the system (7) we prescribe the initial data m(t = 0) = m 0 and ϕ(t = 0) = ϕ 0 in Ω. (2) Additionally, we impose the following physically motivated assumption which we discuss in the next section, 0 The structure of system (7) (posed for Ω = R d ) is derived rigorously in [24] as the hydrodynamic limit of the Kawasaki dynamics for the Blume-Capel model with range of interaction γ −1 , magnetic field h 1 , and chemical potential h 2 , whose Hamiltonian in a finite square where for all r ∈ R d and σ(x) ∈ {−1, 0, +1} is the spin of a particle at position x ∈ V .We discuss the intuition on the discrete spin model below, but also refer the reader to the monograph [31] for more information on the context.As discussed in [24], letting ⃗ u := (m, ϕ), this system can also be written as a gradient flow structure given by where the mobility, and the free energy functional, F, is given by where Also discussed in [24] is that F is a Lyapunov functional for (6) and that Unfortunately, such a variational structure is lost when one makes adjustments to include evaporation.The first equation in system (7) is similar to the non-local Cahn-Hilliard equation well studied in the literature [15,11].In fact, it is easily seen for ϕ ≡ 1 (i.e., a two species system), the model presented here reduces to the non-local Cahn-Hilliard equation with degenerate mobility and logarithmic free energy.However, the coupling of the non-linearity to the second equation, allows for a mechanism to track the solvent ratio in the mixture.Additionally, there is hidden structure in the equations (which we see later in system (9)) that is reminiscent of equations in various applications including the McKean-Vlasov equations studied in swarms [3] and the Fokker-Planck equation studied in mean field games [21] and opinion dynamics [5,4].The main difference between these examples and the equations in system (9), is the non-linearity present in the drift term.While this structure is not surprising to find, as the previously mentioned McKean-Vlasov equations have been used to study interacting particle systems in different applications (e.g., [22]), this particular form of the non-linearity in the drift appears to be quite novel.

Motivating applications and microscopic interpretation
Phase separation models, such as the Cahn-Hilliard equation, have seen a wide variety of applications (see, e.g., [25]).The choice of the particular model ( 7) is motivated by several applications such as the formation of rubber-based zones related to the design of thin adhesive bands [7,26]; or the formation of morphologies in organic solar cells [19].In both of these applications, a mixture consisting of three species of particles (two solutes and one solvent) undergoes phase separation to form the aforementioned morphologies/zones.In both applications, evaporation of the solvent plays a key role as the ratio of solvent to the mixture has been shown to affect the resulting structure of the morphologies (see, e.g., [28] for a discrete simulation and [23] for simulations on system (7)).
Modeling the evaporation process in conjunction with phase separation is usually done at the discrete level with a Monte-Carlo simulation with the solvent particles taking some particular behavior (e.g.,point-wise replacement or replacement at the boundary [6,28,34]).The choice of solvent behavior and spatial dimension changes the interpretation of this evaporation process.Generally for two dimensional simulations, one takes either a 'from the top' or a 'from the side' perspective.The main difficulty is finding the continuum limit of such evaporation models.Therefore, one of the goals of this manuscript is to find a way to adjust system (7), a system derived without an evaporation process, to reasonably simulate evaporation with either perspective.
Discrete lattice simulations of this behavior often assign a 'spin' to the particles denoted as in (4) by σ(x) ∈ {−1, 0, +1}.In this set-up, the solvent particle would be assigned spin 0, while the solute particles would have spin +1 or −1.The concentration of the solution at a site x would then be represented by |σ(x)| (traditionally in the literature as σ(x) 2 ).This is the measure of how much solute (independent of the spin direction) is located at a particular site.
Relating this intuition back to system (7), we let m and ϕ be continuous measures of the average spin (also referred to as magnetization in [24]) and concentration in a given region, respectively.In the context of the discrete model, these terms mean for a collection of sites Therefore at this level, inequalities (3) and ( 8) are just applications of the triangle inequality.Since system (7) is the continuum limit of this stochastic interpretation under the Hamiltonian (4), we maintain this inequality with assumption (3).
Returning back to the continuous setting and system (7), m : QT −→ R is a function describing the net spin of particles in a given region.In other words, for a given Borel set A ⊆ Ω, A m(t, x) dx is the net spin of particles in the set A at time t.As the spin in the Hamiltonian (4) only takes values in {−1, 0, +1}, we aim to show that |m| ≤ 1 almost everywhere in Q T .To capture the nuances of solvent particles, the function ϕ : QT −→ R represents the concentration of solute particles, i.e., A (1 − ϕ(t, x)) dx represents the solvent ratio in the set A at time t.Since ϕ represents the concentration, we will show that the physical inequality 0 ≤ |m| ≤ ϕ ≤ 1 holds almost everywhere in Q T .

A model with evaporation
Keeping all previous assumptions, we now adjust system (1) to include a source term for the concentration component.
Here, F : [0, 1] −→ R is a bounded non-increasing Lipschitz continuous function with F (1) = 0. Since we do not know a priori 0 ≤ ϕ ≤ 1, we extend F to all of R by 0. This reduces the regularity of F , however, the only point of discontinuity is at 0, a fact we exploit later.As F is a source term for the solute concentration function, ϕ, in 2 dimensions, this can represent a 'from the top' evaporation rate.
Our main result is the well-posedness of system (7) and the fact that the qualitative property ( 3) is conserved on Q T .Particularly, we are interested in weak solutions to (7): 2 is a weak solution to (7) if (m(t = 0), ϕ(t = 0)) = (m 0 , ϕ 0 ) and for any t ∈ (0, T ), Here and throughout the paper, we use the subscript ♯ to restrict function spaces to the subset of spatially periodic functions.We aim to prove the following theorem: Theorem 1.1.Let m 0 , ϕ 0 be such that assumption (3) holds.Then, there exists a unique weak solution (m, ϕ) to (7) in the sense of Definition 1.1.Moreover, To this end, assume for a moment that solutions to (7) exist.Then the functions w = ϕ + m and v = ϕ − m would satisfy the system This manipulation has revealed some hidden symmetry in our problem which turns out to be exploitable.Notice now that assumption (3) implies certain bounds on the initial data w 0 , v 0 , specifically, It is easy to see that existence of solutions to system (9) implies the existence of solutions to the original problem (7).Moreover, the linear nature of the transformations between the systems allow us to carry over regularity properties on (w, v) to (m, ϕ).Also, the conservation of inequality (10) for all time will imply the conservation of part of (3), namely, 0 ≤ |m| ≤ ϕ.This combined with an a priori result proven later will yield the conservation of (3) for all t ∈ [0, T ).
Equation ( 9) is written in a peculiar way to foreshadow our plan of attack.We aim to use a Schauder fixed-point argument to prove the existence of solutions by finding a fixed-point on the composite map Such arguments are common in literature, see e.g., [15] for an application to a similar equation to the first in (7) and [9] for an application to systems.

Notation and organization of the paper
Positive constants whose explicit value is unimportant to the story will be denoted arbitrarily by C. If the constant depends on a particular function, f , this will be denoted by C f .We warn that constants who share the same notation, need not share the same value.
We liberally make use of Young's inequality for convolutions: whenever the above norms make sense for the functions in question.
For lucidity, we use a short hand notation for spaces which appear frequently throughout the manuscript.For a given Ω ⊂ R, we define the following sets: The rest of paper is organized in the following way: in Section 2, we show the well-posedness of an auxiliary problem via an iteration scheme; in Section 3 we prove the existence of a fixed-point to equation (9) and complete the proof of Theorem 1.1; in Section 4 we describe a special case of how one would model the evaporation process and provide some numerical simulations.Finally, in Section 5 we conclude with some future work and some discussion on the results.

Analysis of auxiliary problems
Before we can show the well-posedness of system (9) (and therefore system ( 7)), we collect some useful auxiliary and a priori results.First, while a priori boundedness of solutions to (7) is not immediately clear, we can infer that ϕ is bounded given that m is bounded.This result will be useful to us later in showing the map T : m → 1 2 ( w − v) maps bounded functions to bounded functions.Proposition 2.1.Let (m, ϕ) be solutions to (7) in the sense of Definition 1.1.If |m| ≤ 1 almost everywhere in Q T , then it also holds ϕ ≤ 1 almost everywhere in Q T .
Proof.Setting η = [1 − ϕ] − the weak form of ϕ from Definition 1.1, we see that due to the support of F (ϕ), F (ϕ)η = 0.And so, taking the absolute value of the right-hand side of the weak form yields: Using Hölder and Young's inequalities, we arrive at A standard application of the Grönwall inequality (see, e.g., [10,Chapter 7.1]) and using the fact that With similar arguments and motivation, we show the following auxiliary result: Proof.From the weak form of the equation, we have for any η ∈ L 2 (0, T ; Letting η = [u − C] + , we have from the above, By the assumptions on J, we have ∇J * C = 0, we see that ∇J * u = ∇J * (u − C).Therefore, using Hölder and Young's inequalities, A straightforward application of the Grönwall's inequality implies that ∥η∥ 2 2 = 0 for all t ∈ (0, T ).Repeating this argument for the test function ψ = [u + C] − implies the result.

Auxiliary system
In this section, we prove the following system is well-posed: equipped with periodic boundary conditions, where for i = 1, 2, B i : and w(0, x) = w 0 (x), v(0, x) = v 0 (x) are bounded and nonnegative for all x ∈ Ω.Furthermore, we set We understand solutions to equation (12) in the standard weak sense: Inspired by [5,36], we use an iteration scheme method to prove the following theorem Theorem 2.1.There exist a unique non-negative weak solution in the sense of Definition 2.1 to (12) with the following estimates and where ⃗ B := (B 1 , B 2 ).Moreover, Remark.While equations in system (12) are quite similar to the equation studied by authors in [5], they fix a particular function for their convolution and do not consider additional terms in the drift.This allows them certain estimates which are not applicable in our context.

Iteration scheme
We define the iteration scheme as follows: where the initial conditions w(0, x) = w 0 (x) and v(0, x) = v 0 (x) are taken as the initial iteration step.We inherit all of the assumptions from equation (12) and define solutions to this system in a similar manner to Definition 2.1 with the corresponding weak form: We begin with an a priori result: Lemma 2.1.Let {w n } and {v n } be sequences of solutions to (13).Then for each n, and Proof.We prove this statement using mathematical induction.By the assumptions on the data w 0 and v 0 , the claim trivially holds.Now assume w n , v n ≥ 0 with Testing the equations in (14) with ψ = [w n+1 ] − and η = [v n+1 ] − , and summing them together, we have Noticing due to the fact that F is nonegative, we have Therefore, we have the inequality Making use of Hölder's inequality, Young's inequality, and Young's inequality for convolution we have, Taking note that ψ(0, x) ≡ 0 and η(0, x) ≡ 0, the Grönwall inequality implies that ψ ≡ 0 and η ≡ 0 in Q T and, therefore, w n+1 , v n+1 ≥ 0.
Testing equation ( 14) with ψ, η = 1 yields Since w n+1 (0, x) = w 0 (x) and v n+1 (0, x) = v 0 (x) we have and Due to Young's convolution inequality (11), we have that the drift term in ( 13) is in fact uniformly bounded.Indeed, Therefore, inductively applying standard Galerkin arguments found in [10, Section 7.1], we have for each n a unique solution pair (w n+1 , v n+1 ) ∈ X 2 in the sense of ( 14) and with energy estimates similar to those found in Theorem 2.1.Specifically, we have the energy estimates: and for each n = 1, 2, . . . .We now ensure that our scheme is convergent by proving that it is a Cauchy sequence in a suitable space of functions.Notice that since we have shown w n , v n ≥ 0, we may use the Lipschitz continuity of F .Lemma 2.2.There exists a t * > 0 such that the sequence Proof.Let wn := w n − w n−1 and vn := v n − v n−1 .Then from (14) we have that wn+1 and vn+1 satisfy the weak form: The above formulation allows us by testing with (ψ, η) = ( wn+1 , vn+1 ) to arrive at the identities ) Working with the right-hand side of the first equation above we have the first term, , where in the last inequality, we have used Hölder and Young's inequality for convolution.Making use of the weighted Young's inequality and the energy estimate on w n , we can arrive at Looking now to the second term, we have via Young's inequality for convolutions, Hölder's inequality, and Young's product inequality Finally, from the third term on the right-hand side we see by Lipschitz continuity of F (denoting the Lipschitz constant by L F ), , where in the last step we again make use of Hölder's inequality and Young's inequality.
Repeating these calculations for vn+1 and summing the equations in (16), we arrive at the estimates Choosing t * such that C ⃗ B,J,F,w0,v0 t * 1 − C ⃗ B,J,F,w0,v0 t * < 1 and integrating (17) from 0 to t < t * we can arrive at the estimate Taking the supremum of both sides over [0, t * ), we obtain: where we have used the inequality Manipulation of the above inequality yields which implies the sequence is Cauchy in L ∞ (0, t * ; L 2 (Ω)).Moreover, using inequality (18), we see this implies the sequence is also Cauchy in L 2 (0, t * ; L 2 (Ω)).
Remark.The uniqueness of solutions to (12) follows from similar arguments presented in the above proof.Indeed, letting w = w 1 − w 2 and v = v 1 − v 2 where (w 1 , v 1 ) and (w 2 , v 2 ) are solutions to (12), we see that w and v satisfy similar equations to (15), namely Following the calculations which lead to equation (17), we arrive at the differential inequality A standard application of Grönwall's inequality implies ∥ w∥ 2 2 + ∥v∥ 2 2 = 0 for all t ∈ [0, T ].Remark.By continuity and the above energy estimates, we can extend the solution to any finite T using standard methods (see, e.g., [1,Section 7]).
Summing up, we have obtained that the constructed sequence (w n , v n ) is in fact a Cauchy sequence in L ∞ (0, T ; L 2 (Ω)) 2 (or L 2 (Q T ) 2 ) and, from the energy estimates, is weakly convergent (along a subsequence) in L 2 (0, T ; Using these results, we can pass to the limit in the weak formulation ( 14) (following the calculations used to derive equation ( 17)) as n −→ ∞ to arrive at the weak formulation of Definition 2.1.Thus, the limit is a solution to (12).

The fixed-point argument
We are now ready to prove the well-posedness of equation ( 9) which for the reader's convenience we rewrite here: Taking into account assumption (10), solutions to (9) are understood as in Definition 2.1 with B 1 := β(m − 1) and B 2 := β(m + 1).We aim to show the existence of solutions to a fixed-point to the map where ) is the weak solution to (12) with B 1 = β( m − 1) and B 2 = β( m + 1), and T 2 ( w, v) = 1 2 ( w − v).Now, using the work done in the previous sections, we have the following lemmas: Lemma 3.1.The map T = T 2 • T 1 is well-defined over the set Z and T : Z −→ Z.
Proof.We study separately the maps T 1 and T 2 .

The map T 1 :
Fix m ∈ Z.By Theorem 2.1, we have the existence of a unique pair 2 such that 0 ≤ w, v. Thus the map T 1 is well-defined.
The map T 2 : The map T 2 is trivially well-defined.What remains to be shown is that T 2 maps into Z.To see this, define with . Thus, by Proposition 2.2, we have that |m| ≤ 1.
Remark 3.1.With similar arguments to those found in Proposition 2.1, it can be easily shown the function ϕ := 1 2 ( w + v), which satisfies the equation is nonnegative and bounded by 1 (i.e.,0 ≤ ϕ ≤ 1) since m, m ∈ Z and 0 ≤ w, v. Together, these bounds allow us to conclude 0 ≤ w, v ≤ 2 almost everywhere in Q T .
Lemma 3.2.The map Proof.Let mn , m ∈ Z be such that mn −→ m in L 2 (Q T ).Define also ( wn , vn ) := T 1 (m n ) and ( w, v) := T 1 (m).We want to show ( wn , vn ) −→ ( w, v) in L 2 (Q T ).To this end, let m := mn − m, w := wn − w, and v := vn − v for fixed n.Following similar manipulations which gave rise to equation (15), we have that w and v satisfy Testing with (ψ, η) = (w, v) we have the following equations: Making use of Remark 3.1, the first two terms on the right-hand side of both equations can be handled as in the proof of Lemma 2.2 to arrive at the bounds: Similarly, the evaporation terms follow the same arguments, Meanwhile, the new third term yields the bound Summing up the two inequalities and maximizing constants leaves us with Integrating the above equation from 0 to t gives us the integral inequality Grönwall's inequality and integration provides us with the estimate: Continuity of T 1 and therefore T follows.Proof.Summing up the results from before, we have the set Z is a non-empty, closed, and convex subset of L 2 (Q T ), the map T : Z −→ Z is well-defined and continuous (Lemmas 3.1 and 3.2), and T (Z) ⊆ W which is compactly embedded in L 2 (Q T ) via Lions-Aubin's Lemma.Therefore, we can employ Schauder's Fixed-Point Theorem to conclude there exists a fixed-point in W ∩ Z.
Summing up, we have shown the existence of bounded solutions (w, v) to problem (9).By manipulation, we can see that the functions m = 1 2 (w − v) and ϕ = 1 2 (w + v) are solutions to equation (7).Moreover, we have shown that 0 ≤ w, v ≤ 2 which implies that |m| ≤ ϕ.
Since in the fixed-point argument we show that m ∈ Z, we can conclude via Proposition 2.1 that ϕ ≤ 1.Thus we have shown the existence of bounded solutions to equation (7) in the sense of Definition 1.1.A proof of the uniqueness of bounded solutions to system (1) can be found in [24,Lemma 5.4].However, one can also prove uniqueness in L 2 for system (7) with similar arguments to Lemma 3.2 above.Indeed, we have Lemma 3.3.The solution pair (w, v) given by Definition 2.1 is unique.
Proof.Letting w = w 1 − w 2 and v = v 1 − v 2 be the difference of two solutions of (7), following the calculations in Lemma 3.2, and using that m = 1 2 (w − v), we have that inequality (26)  Grönwall's inequality then proves the result.

A simple example of 'from the top' evaporation
Looking at a thin film, modeling the evaporation of the solvent phase from a mixture of multiple phases interacting among themselves as well as with the solvent, is, in general, quite complicated.Consequently, accurate quantitative predictions are usually out of reach.However, in an experimentally controlled setup with not too many interacting phases and where the effect of temperature gradients can be neglected, valuable qualitative insights can be reached.For instance, in the context of applications to organic solar cells, the solvent evaporates into vapor and moves up and leaves the film.This process leaves us with two perspectives on the actual geometry of the thin film: from 'the side' and from 'the top'.For continuum models, the side perspective of the evaporation process can be modeled as a moving interface separating the solvent vapor and liquid phases as done, for instance, in [8].This process tends to result in the formation of solvent 'lanes' and can be observed both in continuum [8] models as well as in discrete lattice models [34].On the other hand, an observer of the top view of the film would see a bulk evaporation which is usually modeled as right-hand side source term of a balance law for solvent posed in two dimensions and not as a moving-boundary condition imposed an a moving surface within the thin film.Such situation is conceptually simpler and is easier to handle from modeling, mathematical analysis, and simulation points of view.In this section, we use a simple finite volume method (previously tested on model (1) in [23]) to demonstrate 'from the top' evaporation and the effect it has on the morphologies.We leave the 'from the side' and the 3 dimensional combination of the two perspectives for future work.

Simulation results
To fix ideas, we take into account a linear evaporation model, i.e., F : R → R defined by F (r) = α(1−r) for all r ∈ R with given α > 0. We refer the reader to section 2.4.3.3 in [27] for a motivation of this particular structure that describes in a very simplified way the liquid-gas transition.More information on how to include the evaporation mechanism in models capturing phase separation in thin films can be found e.g. in [33].
We use a finite volume scheme similar to that presented in [23] with slight modifications to include the production term by evaporation.Namely, for fixed mesh sizes ∆t, ∆x, ∆y > 0 we discretize the domain Ω by cells Λ i,j := [x i − 1 2 h, x i + 1 2 h) × [y j − 1 2 h, y j + 1 2 h) for i, j = 1, 2, . . ., N .To account for the periodic boundary conditions, we periodically extend functions defined on the nodes, i.e., f i,j = f i±N,j = f i,j±N .The initial pair (m 0 , ϕ 0 ) is then approximated by m 0 (x, y) dx dy and ϕ 0 i,j := 1 |Λ i,j | Λi,j ϕ 0 (x, y) dx dy and define the fully explicit scheme:

and Jk
x,i,j (respectively, Jk y,i,j ) denotes the approximation of ∂ x J * m k (x i , y j ) (respectively, ∂ y J * m k (x i , y j )) using a fast Fourier transform method similar to [35].We remark that similar estimates to those found in Section 2 may be repurposed to discuss convergence of the finite volume scheme.However, we leave this to be rigorously shown in a future manuscript.A numerical demonstration of convergence for the scheme with F (ϕ) ≡ 0 can be found in [23].
As in [23], we observe morphology formation during the simulation.However, whereas in the previous results the shape-type (e.g.,'ball-like' or 'bi-continuous') of morphologies stays consistent throughout the simulation, we observe a change in shape-type while there is an evaporation component.This phenomenon is demonstrated in Figure 1 where 'ball-like' structures can be seen at early times (t = 1) and they slowly grow as solvent evaporates into long continuous structures at later times.A similar situation pointing out a competition between phase separation and evaporation is reported in [32].Here the shapes of the morphologies formed in a partially miscible mixture change during the evaporation process, while one notes that the spinodal instability occurs when the evaporation is fast.In Figure 2, we present multiple measurements of the simulation.First, we present the ratio of solvent during the simulation over time given by ∥1 − ϕ∥ 1 .As expected, we observe a decrease in solvent ratio over time.Next, we plot the ∥ • ∥ 1 norms of both m and ϕ (scaled by the norm of the initial condition) over time.The observed behavior is an increase in total volume of the solute.This is reminiscent of the Monte Carlo simulation of somewhat related microscopic models (see [28,6]), where the evaporating solvent was replaced by solute particles.We note that the sudden drop in the L 1 norm of m is due to the initial mixing of the solution (ternary mixture).In other words, this drop happens during the time period where morphologies are not yet well formed.In the aforementioned cited papers, the evaporating solvent was replaced by solute in such a way that the ratio of +1 and −1 particles was kept constant.We measure this effect in the third plot of Figure 2 where we plot the evolution of the ratio of ∥m + ∥ 1 and ∥m − ∥ 1 .Here, we observe an approximate 4% variance of this ratio due to the numerical scheme.However, it appears to self correct over time.

Discussion and outlook
Summarizing our work here, we have shown the existence of a unique bounded solution to system (7) with standard parabolic regularity for some finite time.Moreover, the solutions to this system conserve the physical inequality (8).We have included an example of how to adjust system (7) to include a 'from the top' evaporation process and included simulations of a simple example of this process.While the particular example of evaporation is simple, clear transitions between the patterns studied in [23] can be observed in Figure 1.
The study of such type of phase separation models is rich with questions and it remains to be seen how the analysis of similar models such as the Cahn-Hilliard and Allen-Cahn equations compare to system (7).For instance, the study of the sharp interface via formal asymptotic expansion as done in [15,30,18,2] could lead to more insight on how the morphologies grow in time.One could add further complications to this study by applying the asymptotic expansion to system (7) including the evaporation process.Another question relating to the analysis of the model is with respect to the strict separation property studied in [11,13,12].Numerical simulations suggest that ||m| − ϕ| −→ 0 almost everywhere as t −→ ∞ for initial conditions satisfying assumption (3).This property has the physical meaning that the mixture asymptotically becomes well separated and has yet to be rigorously shown.
Our main interest is to use the morphologies generated by system (7) as a type of porous domain to model and simulate charge transport under some accepted physical assumptions [29,20].Such simulations can allow a quantitative study of the effectiveness of the patterns of observed morphologies.With this in mind, more simulations of system (7) in three dimensions taking both perspectives of evaporation into account are crucial to best representing the physical system.Numerical schemes taking into account the gradient flow structure of the model presented in [24] similar to those studied in [17,16] seem readily adaptable to system (7) without evaporation.As future work, we would like to study ways to include the evaporation process inside of these schemes.
with periodic boundary conditions and where u(t = 0) = u 0 with |u 0 | ≤ C almost everywhere in Ω.Then, |u| ≤ C almost everywhere in Q T .

Figure 1 :
Figure 1: Simulation of equation (7) with α = 0.1, β = 10, and initial solvent ratio of 80%.Regions where m is positive are colored blue, negative regions are colored yellow, and regions m is near zero are colored red.Similarly, regions where ϕ is near one are colored blue whereas regions near zero are colored red.

Figure 2 :
Figure 2: Solvent ratio, L 1 norms, and solute ratio of the simulation shown in Figure 1 plotted over time.The horizontal axis here represents the time step in the numerical scheme in thousands.