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

Rainey<EMAIL_ADDRESS>Emilio N. M. Cirillo  Adrian Muntean

###### Abstract

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 ‘from
the top’ 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.

Keywords: Phase separation, coupled non-local parabolic system, ternary
mixture, evaporation, well-posedness, weak solutions, numerical simulation AMS
Subject classifications: 35K55, 35Q70, 65N08

##  1 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–Cahn–type 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.

###  1.1 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

$\left\\{\begin{split}\partial_{t}m&=\nabla\cdot\left[\nabla
m-2\beta(\phi-m^{2})(\nabla
J*m)\right],\quad(t,x)\in(0,T)\times\Omega=:Q_{T},\\\
\partial_{t}\phi&=\nabla\cdot\left[\nabla\phi-2\beta m(1-\phi)(\nabla
J*m)\right].\end{split}\right.$ (1)

Here, $t$ and $x$ represent the time and space variable, respectively, $m$
represents the average spin density (also called magnetization), and $\phi$
represents the solute volume concentration (we expand upon the meaning of
these terms in a later section). Additionally, $\Omega\,\subset\mathbb{R}^{d}$
is a cube with spatially periodic boundary conditions (i.e., homeomorphic to
the unit torus), $T$ is some positive time, $\beta>0$ is a constant
representing the inverse temperature, $J\in C_{+}^{2}(\bar{\Omega})$ is a
symmetric compactly supported potential such that $\int_{\Omega}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\in W^{1,\infty}(\mathbb{R}^{d})$
and small adjustments can be made for to include $J\in
W^{1,2}(\mathbb{R}^{d})$ and possibly $J\in W^{1,1}(\mathbb{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}\mbox{ and }\phi(t=0)=\phi_{0}\mbox{ in }\bar{\Omega}.$ (2)

Additionally, we impose the following physically motivated assumption which we
discuss in the next section,

$0\leq|m_{0}|\leq\phi_{0}\leq 1\text{ almost everywhere in }\bar{\Omega}.$ (3)

The structure of system (7) (posed for $\Omega=\mathbb{R}^{d}$) is derived
rigorously in [24] as the hydrodynamic limit of the Kawasaki dynamics for the
Blume–Capel model with range of interaction $\gamma^{-1}$, magnetic field
$h_{1}$, and chemical potential $h_{2}$, whose Hamiltonian in a finite square
$V\subset\mathbb{Z}^{d}$ is

$H_{\gamma}(\sigma)=\frac{1}{2}\sum_{\begin{subarray}{c}x\neq x^{\prime}\in
V\end{subarray}}J_{\gamma}(x-x^{\prime})[\sigma(x)-\sigma(x^{\prime})]^{2}-\sum_{x\in
V}h_{1}\sigma(x)-\sum_{x\in V}h_{2}\sigma^{2}(x),$ (4)

where $J_{\gamma}:\mathbb{R}^{d}\to\mathbb{R}$ is a Kac potential function,
i.e.,

$J_{\gamma}(r)=\gamma^{d}J(\gamma r)$ (5)

for all $r\in\mathbb{R}^{d}$ and $\sigma(x)\in\\{-1,0,+1\\}$ is the spin of a
particle at position $x\in 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 $\vec{u}:=(m,\phi)$, this system can also be
written as a gradient flow structure given by

$\partial_{t}\vec{u}=\nabla\cdot\left(M\nabla\frac{\delta\mathcal{F}}{\delta\vec{u}}\right),$
(6)

where the mobility,

$M:=\beta(1-\phi)\begin{bmatrix}\phi+\frac{\phi^{2}-m^{2}}{1-\phi}&m\\\
m&\phi\end{bmatrix}$

and the free energy functional, $\mathcal{F}$, is given by

$\mathcal{F}(\vec{u})=\int_{\Omega}f(\vec{u})\mathop{}\\!\mathrm{d}{x}+\frac{1}{2}\int_{\Omega}\int_{\Omega}J(x-x^{\prime})[m(x)-m(x^{\prime})]^{2}\mathop{}\\!\mathrm{d}{x^{\prime}}\mathop{}\\!\mathrm{d}{x},$

where

$f(\vec{u}):=\phi-m^{2}+\beta^{-1}\left[\frac{1}{2}(\phi+m)\log(\phi+m)+\frac{1}{2}(\phi-m)\log(\phi-m)+(1-\phi)\log(1-\phi)-\phi\log(2)\right].$

Also discussed in [24] is that $\mathcal{F}$ is a Lyapunov functional for (6)
and that

$\frac{\mathop{}\\!\mathrm{d}{}}{\mathop{}\\!\mathrm{d}{t}}\mathcal{F}=-\sum_{i=1}^{d}\sum_{\alpha,\beta=1}^{2}\int_{\Omega}\partial_{i}\frac{\delta\mathcal{F}}{\delta\vec{u}_{\alpha}}\,M_{\alpha,\beta}\,\partial_{i}\frac{\delta\mathcal{F}}{\delta\vec{u}_{\beta}}.$

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 $\phi\equiv 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.

###  1.2 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 $\sigma(x)\in\\{-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 $|\sigma(x)|$ (traditionally in the literature as
$\sigma(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 $\phi$ 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
$A:=\\{x_{i}\\}_{i=1}^{N}$

$m_{A}=\frac{1}{N}\sum_{i=1}^{N}\sigma(x_{i})\text{ and
}\phi_{A}=\frac{1}{N}\sum_{i=1}^{N}\sigma(x_{i})^{2}=\frac{1}{N}\sum_{i=1}^{N}|\sigma(x_{i})|.$

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:\bar{Q}_{T}\longrightarrow\mathbb{R}$ is a function describing the net spin
of particles in a given region. In other words, for a given Borel set
$A\subseteq\Omega$, $\int_{A}m(t,x)\mathop{}\\!\mathrm{d}x$ 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|\leq 1$ almost
everywhere in $Q_{T}$. To capture the nuances of solvent particles, the
function $\phi:\bar{Q}_{T}\longrightarrow\mathbb{R}$ represents the
concentration of solute particles,
i.e.,$\int_{A}(1-\phi(t,x))\mathop{}\\!\mathrm{d}x$ represents the solvent
ratio in the set $A$ at time $t$. Since $\phi$ represents the concentration,
we will show that the physical inequality $0\leq|m|\leq\phi\leq 1$ holds
almost everywhere in $Q_{T}$.

###  1.3 A model with evaporation

Keeping all previous assumptions, we now adjust system (1) to include a source
term for the concentration component.

$\left\\{\begin{split}\partial_{t}m&=\nabla\cdot\left[\nabla
m-2\beta(\phi-m^{2})(\nabla
J*m)\right],\quad(t,x)\in(0,T)\times\Omega=:Q_{T},\\\
\partial_{t}\phi&=\nabla\cdot\left[\nabla\phi-2\beta m(1-\phi)(\nabla
J*m)\right]+F(\phi).\end{split}\right.$ (7)

Here, $F:[0,1]\longrightarrow\mathbb{R}$ is a bounded non-increasing Lipschitz
continuous function with $F(1)=0$. Since we do not know a priori
$0\leq\phi\leq 1$, we extend $F$ to all of $\mathbb{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,
$\phi$, 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):

######  Definition 1.1.

We say a pair $(m,\phi)\in\left(L^{2}(0,T;H^{1}_{\sharp}(\Omega))\cap
L^{\infty}(0,T;L^{2}(\Omega))\right)^{2}$ with
$(\partial_{t}m,\partial_{t}\phi)\in\left(L^{2}(0,T;H^{-1}_{\sharp}(\Omega))\right)^{2}$
is a weak solution to (7) if $(m(t=0),\phi(t=0))=(m_{0},\phi_{0})$ and for any
$t\in(0,T)$,

$\langle\partial_{t}m,\psi\rangle+\int_{\Omega}\nabla
m\cdot\nabla\psi-2\beta(\phi-m^{2})(\nabla
J*m)\cdot\nabla\psi\mathop{}\\!\mathrm{d}x=0,$

and

$\langle\partial_{t}\phi,\eta\rangle+\int_{\Omega}\nabla\phi\cdot\nabla\eta-2\beta
m(1-\phi)(\nabla
J*m)\cdot\nabla\eta\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\phi)\eta\mathop{}\\!\mathrm{d}{x},$

for all $(\psi,\eta)\in L^{2}(0,T;H^{1}_{\sharp}(\Omega))\times
L^{2}(0,T;H^{1}_{\sharp}(\Omega))$.

Here and throughout the paper, we use the subscript $\sharp$ to restrict
function spaces to the subset of spatially periodic functions. We aim to prove
the following theorem:

######  Theorem 1.1.

Let $m_{0},\phi_{0}$ be such that assumption (3) holds. Then, there exists a
unique weak solution $(m,\phi)$ to (7) in the sense of Definition 1.1.
Moreover,

$0\leq|m|\leq\phi\leq 1\text{ almost everywhere in }Q_{T}.$ (8)

To this end, assume for a moment that solutions to (7) exist. Then the
functions $w=\phi+m$ and $v=\phi-m$ would satisfy the system

$\left\\{\begin{split}\partial_{t}w&=\Delta w+\text{div}[\beta(m-1)(\nabla
J*(w-v))w]+F(\tfrac{1}{2}(w+v)),\\\ \partial_{t}v&=\Delta
v+\text{div}[\beta(m+1)(\nabla J*(w-v))v]+F(\tfrac{1}{2}(w+v)),\\\
m&=\frac{1}{2}(w-v),\end{split}\right.\quad(t,x)\in Q_{T}.$ (9)

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,

$0\leq w_{0},v_{0}\leq 2.$ (10)

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,\phi)$. Also, the conservation of
inequality (10) for all time will imply the conservation of part of (3),
namely, $0\leq|m|\leq\phi$. This combined with an a priori  result proven
later will yield the conservation of (3) for all $t\in[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

$\mathcal{T}:\bar{m}\xmapsto{\mathcal{T}_{1}}(\bar{w},\bar{v})\xmapsto{\mathcal{T}_{2}}\frac{1}{2}(\bar{w}-\bar{v}).$

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.

###  1.4 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:

$\|f*g\|_{r}\leq\|f\|_{p}\|g\|_{q},\quad\text{where}\quad\frac{1}{p}+\frac{1}{q}=1+\frac{1}{r},$
(11)

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 $\Omega\subset\mathbb{R^{d}}$, we
define the following sets:

$\mathcal{X}:=L^{2}(0,T;H^{1}_{\sharp}(\Omega))\cap
L^{\infty}(0,T;L^{2}(\Omega));$ $\mathcal{W}:=\left\\{u\in
L^{2}(0,T;H^{1}_{\sharp}(\Omega)):\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}{t}}u\in
L^{2}(0,T;H^{-1}_{\sharp}(\Omega))\right\\}.$

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.

##  2 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 $\phi$ is bounded given that $m$ is bounded. This result will be useful
to us later in showing the map
$\mathcal{T}:\bar{m}\mapsto\frac{1}{2}(\bar{w}-\bar{v})$ maps bounded
functions to bounded functions.

######  Proposition 2.1.

Let $(m,\phi)$ be solutions to (7) in the sense of Definition 1.1. If $|m|\leq
1$ almost everywhere in $Q_{T}$, then it also holds $\phi\leq 1$ almost
everywhere in $Q_{T}$.

###### Proof.

Setting $\eta=[1-\phi]^{-}$ the weak form of $\phi$ from Definition 1.1, we
see that due to the support of $F(\phi)$, $F(\phi)\eta=0$. And so, taking the
absolute value of the right-hand side of the weak form yields:

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\eta\|_{2}^{2}+\|\nabla\eta\|_{2}^{2}\leq\int_{\Omega}|(2\beta
m\eta\nabla J*m)\cdot\nabla\eta|\mathop{}\\!\mathrm{d}x.$

Using Hölder and Young’s inequalities, we arrive at

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\eta\|_{2}^{2}+\|\nabla\eta\|_{2}^{2}\leq
2\beta^{2}C_{J}\|\eta\|_{2}^{2}+\frac{1}{2}\|\nabla\eta\|_{2}^{2}.$

A standard application of the Grönwall inequality (see, e.g., [10, Chapter
7.1]) and using the fact that $\eta(t=0)\equiv 0$ implies that
$\eta=[1-\phi]^{-}=0$ almost everywhere in $Q_{T}$. ∎

With similar arguments and motivation, we show the following auxiliary result:

######  Proposition 2.2.

Assume for some $B:(0,T)\times\Omega\longrightarrow\mathbb{R}$ with
$\|B\|_{L^{\infty}(0,T;L^{2}(\Omega))}\leq C_{B}$ and $J\in
C_{+}^{2}(\Omega)$, we have a weak solution $u\in
L^{2}(0,T;H^{1}_{\sharp}(\Omega))$ to the equation

$\partial_{t}u=\Delta u+\text{div}(B\nabla J*u),\quad(t,x)\in Q_{T}$

equipped with periodic boundary conditions and where $u(t=0)=u_{0}$ with
$|u_{0}|\leq C$ almost everywhere in $\bar{\Omega}$. Then, $|u|\leq C$ almost
everywhere in $Q_{T}$.

###### Proof.

From the weak form of the equation, we have for any $\eta\in
L^{2}(0,T;H^{1}_{\sharp}(\Omega))$,

$\langle\partial_{t}u,\eta\rangle+\int_{\Omega}\nabla
u\cdot\nabla\eta+(B\nabla J*u)\cdot\nabla\eta\mathop{}\\!\mathrm{d}x=0.$

Letting $\eta=[u-C]^{+}$, we have from the above,

$\frac{d}{dt}\|\eta\|_{2}^{2}+\|\nabla\eta\|_{2}^{2}=-\int_{\Omega}(B\nabla
J*u)\cdot\nabla\eta\mathop{}\\!\mathrm{d}x.$

By the assumptions on $J$, we have $\nabla J*C=0$, we see that $\nabla
J*u=\nabla J*(u-C)$. Therefore, using Hölder and Young’s inequalities,

$\displaystyle\frac{d}{dt}\|\eta\|_{2}^{2}+\|\nabla\eta\|_{2}^{2}$
$\displaystyle\leq\|\nabla J*\eta\|_{\infty}\,\|B\|_{2}\,\|\nabla\eta\|_{2}$
$\displaystyle\leq\|\nabla
J\|_{2}\,\|\eta\|_{2}\,\|B\|_{2}\,\|\nabla\eta\|_{2}$
$\displaystyle\leq\frac{1}{2}C_{J}^{2}\,C_{B}^{2}\,\|\eta\|_{2}^{2}+\frac{1}{2}\|\nabla\eta\|_{2}^{2}.$

A straightforward application of the Grönwall’s inequality implies that
$\|\eta\|_{2}^{2}=0$ for all $t\in(0,T)$. Repeating this argument for the test
function $\psi=[u+C]^{-}$ implies the result.

∎

###  2.1 Auxiliary system

In this section, we prove the following system is well-posed:

$\left\\{\begin{split}\partial_{t}w&=\Delta w+\text{div}[B_{1}(\nabla
J*(w-v))w]+F(\tfrac{1}{2}(w+v))\\\ \partial_{t}v&=\Delta
v+\text{div}[B_{2}(\nabla J*(w-v))v]+F(\tfrac{1}{2}(w+v))\\\
\end{split}\right.,\quad(t,x)\in Q_{T},$ (12)

equipped with periodic boundary conditions, where for $i=1,2$,
$B_{i}:Q_{T}\longrightarrow\mathbb{R}\in L^{\infty}(Q_{T})$ and
$w(0,x)=w_{0}(x)$, $v(0,x)=v_{0}(x)$ are bounded and nonnegative for all
$x\in\bar{\Omega}$. Furthermore, we set
$\int_{\Omega}w_{0}\mathop{}\\!\mathrm{d}x=K_{w}$,
$\int_{\Omega}v_{0}\mathop{}\\!\mathrm{d}x=K_{v}$ for some $K_{w},K_{v}\geq
0$. We understand solutions to equation (12) in the standard weak sense:

######  Definition 2.1.

We say a pair $(w,v)\in\mathcal{X}^{2}$ with $(\partial_{t}w,\partial_{t}v)\in
L^{2}(0,T;H^{-1}_{\sharp}(\Omega))^{2}$ is a solution to (12) if
$(w(t=0),v(t=0))=(w_{0},v_{0})$ and for every $t\in(0,T)$,

$\langle\partial_{t}w,\psi\rangle+\int_{\Omega}\nabla
w\cdot\nabla\psi+wB_{1}(\nabla
J*(w-v))\cdot\nabla\psi\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\tfrac{1}{2}(w+v))\psi\mathop{}\\!\mathrm{d}{x},$

and

$\langle\partial_{t}v,\eta\rangle+\int_{\Omega}\nabla
v\cdot\nabla\eta+vB_{2}(\nabla
J*(w-v))\cdot\nabla\eta\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\tfrac{1}{2}(w+v))\eta\mathop{}\\!\mathrm{d}{x},$

for all $(\psi,\eta)\in L^{2}(0,T;H^{1}_{\sharp}(\Omega))^{2}$.

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

$\|w\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|w\|_{L^{2}(0,T;H^{1}_{\sharp}(\Omega))}+\|\partial_{t}w\|_{L^{2}(0,T;H^{-1}_{\sharp}(\Omega))}\leq
C_{\vec{B},J}(\|w_{0}\|_{2}+\|v_{0}\|_{2}+\|F\|_{\infty}\sqrt{T}),$

and

$\|v\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|v\|_{L^{2}(0,T;H^{1}_{\sharp}(\Omega))}+\|\partial_{t}v\|_{L^{2}(0,T;H^{-1}_{\sharp}(\Omega))}\leq
C_{\vec{B},J}(\|w_{0}\|_{2}+\|v_{0}\|_{2}+\|F\|_{\infty}\sqrt{T}),$

where $\vec{B}:=(B_{1},B_{2})$. Moreover,

$\int_{\Omega}w\mathop{}\\!\mathrm{d}x\leq K_{w}+C_{F}t\text{ and
}\int_{\Omega}v\mathop{}\\!\mathrm{d}x\leq K_{v}+C_{F}t\text{ for all
}t\in[0,T).$

######  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.

###  2.2 Iteration scheme

We define the iteration scheme as follows:

$\left\\{\begin{split}\partial_{t}w_{n+1}&=\Delta
w_{n+1}+\text{div}[B_{1}w_{n+1}\,\nabla
J*(w_{n}-v_{n})]+F(\tfrac{1}{2}(w_{n}+v_{n})),\\\ \partial_{t}v_{n+1}&=\Delta
v_{n+1}+\text{div}[B_{2}v_{n+1}\,\nabla
J*(w_{n}-v_{n})]+F(\tfrac{1}{2}(w_{n}+v_{n})),\\\
\end{split}\right.\quad(t,x)\in Q_{T},$ (13)

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:

$\left\\{\begin{split}\langle\partial_{t}w_{n+1},\psi\rangle+\int_{\Omega}\nabla
w_{n+1}\cdot\nabla\psi+w_{n+1}B_{1}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\psi\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\tfrac{1}{2}(w_{n}+v_{n}))\psi\mathop{}\\!\mathrm{d}{x},\\\
\langle\partial_{t}v_{n+1},\eta\rangle+\int_{\Omega}\nabla
v_{n+1}\cdot\nabla\eta+v_{n+1}B_{2}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\eta\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\tfrac{1}{2}(w_{n}+v_{n}))\eta\mathop{}\\!\mathrm{d}{x}.\end{split}\right.$
(14)

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$,

$w_{n},v_{n}\geq 0,\int_{\Omega}w_{n}\mathop{}\\!\mathrm{d}x\leq
K_{w}+\|F\|_{\infty}|\Omega|t,\text{ and
}\int_{\Omega}v_{n}\mathop{}\\!\mathrm{d}x\leq K_{v}+\|F\|_{\infty}|\Omega|t.$

###### 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}\geq 0$ with

$\int_{\Omega}w_{n}\mathop{}\\!\mathrm{d}x\leq
K_{w}+\|F\|_{\infty}|\Omega|t\text{ and
}\int_{\Omega}v_{n}\mathop{}\\!\mathrm{d}x\leq K_{v}+\|F\|_{\infty}|\Omega|t.$

Testing the equations in (14) with $\psi=[w_{n+1}]^{-}$ and
$\eta=[v_{n+1}]^{-}$, and summing them together, we have

$\displaystyle\frac{1}{2}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\psi\|_{2}^{2}+\|\nabla\psi\|_{2}^{2}+\frac{1}{2}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\eta\|_{2}^{2}+\|\nabla\eta\|_{2}^{2}$
$\displaystyle=-\int_{\Omega}B_{1}\psi\nabla
J*(w_{n}-v_{n})\cdot\nabla\psi\mathop{}\\!\mathrm{d}x-\int_{\Omega}B_{2}\eta\nabla
J*(w_{n}-v_{n})\cdot\nabla\eta\mathop{}\\!\mathrm{d}x$
$\displaystyle\quad+\int_{\Omega}F(\tfrac{1}{2}(w_{n}+v_{n}))(\psi+\eta)\mathop{}\\!\mathrm{d}{x}.$

Noticing due to the fact that $F$ is nonegative, we have

$\int_{\Omega}F(\tfrac{1}{2}(w_{n+1}+v_{n+1}))(\psi+\eta)\mathop{}\\!\mathrm{d}{x}\leq\int_{\Omega}F(\tfrac{1}{2}(w_{n}+v_{n}))[w_{n+1}+v_{n+1}]^{-}\mathop{}\\!\mathrm{d}{x}\leq
0.$

Therefore, we have the inequality

$\frac{1}{2}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\psi\|_{2}^{2}+\|\nabla\psi\|_{2}^{2}+\frac{1}{2}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\eta\|_{2}^{2}+\|\nabla\eta\|_{2}^{2}\leq-\int_{\Omega}B_{1}\psi\nabla
J*(w_{n}-v_{n})\cdot\nabla\psi\mathop{}\\!\mathrm{d}x-\int_{\Omega}B_{2}\eta\nabla
J*(w_{n}-v_{n})\cdot\nabla\eta\mathop{}\\!\mathrm{d}x.$

Making use of Hölder’s inequality, Young’s inequality, and Young’s inequality
for convolution we have,

$\displaystyle\frac{1}{2}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\psi\|_{2}^{2}+\frac{1}{2}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\eta\|_{2}^{2}$
$\displaystyle\leq\frac{1}{2}\|B_{1}\nabla
J*(w_{n}-v_{n})\|_{\infty}^{2}\|\psi\|_{2}^{2}+\frac{1}{2}\|B_{2}\nabla
J*(w_{n}-v_{n})\|_{\infty}^{2}\|\eta\|_{2}^{2}$
$\displaystyle\leq\frac{1}{2}(C_{B}C_{J}(K_{w}+K_{v}+2\|F\|_{\infty}|\Omega|T))^{2}(\|\psi\|_{2}^{2}+\|\eta\|_{2}^{2}).$

Taking note that $\psi(0,x)\equiv 0$ and $\eta(0,x)\equiv 0$, the Grönwall
inequality implies that $\psi\equiv 0$ and $\eta\equiv 0$ in $Q_{T}$ and,
therefore, $w_{n+1},v_{n+1}\geq 0.$

Testing equation (14) with $\psi,\eta=1$ yields

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\int_{\Omega}w_{n+1}\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\tfrac{1}{2}(w_{n}+v_{n}))\mathop{}\\!\mathrm{d}{x}\text{
and
}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\int_{\Omega}v_{n+1}\mathop{}\\!\mathrm{d}x=\int_{\Omega}F(\tfrac{1}{2}(w_{n}+v_{n}))\mathop{}\\!\mathrm{d}{x}.$

Since $w_{n+1}(0,x)=w_{0}(x)$ and $v_{n+1}(0,x)=v_{0}(x)$ we have

$\int_{\Omega}w_{n+1}\mathop{}\\!\mathrm{d}x=\int_{\Omega}w_{0}\mathop{}\\!\mathrm{d}x+\int_{0}^{t}\int_{\Omega}F(\tfrac{1}{2}(w_{n}(s,\cdot)+v_{n}(s,\cdot)))\mathop{}\\!\mathrm{d}{x}\mathop{}\\!\mathrm{d}{t}\leq
K_{w}+\|F\|_{\infty}|\Omega|t$

and

$\int_{\Omega}v_{n+1}\mathop{}\\!\mathrm{d}x=\int_{\Omega}v_{0}\mathop{}\\!\mathrm{d}x+\int_{0}^{t}\int_{\Omega}F(\tfrac{1}{2}(w_{n}(s,\cdot)+v_{n}(s,\cdot)))\mathop{}\\!\mathrm{d}{x}\mathop{}\\!\mathrm{d}{t}\leq
K_{v}+\|F\|_{\infty}|\Omega|t$

for all $t\in(0,T)$. ∎

Due to Young’s convolution inequality (11), we have that the drift term in
(13) is in fact uniformly bounded. Indeed,

$\|B_{i}\nabla
J*(w_{n}-v_{n})\|_{L^{\infty}(Q_{T})}\leq\|B_{i}\|_{L^{\infty}(Q_{T})}\|\nabla
J\|_{\infty}(K_{w}+K_{v}+2|\Omega|\,T\,\|F\|_{\infty}).$

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})\in\mathcal{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:

$\|w_{n}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|w_{n}\|_{L^{2}(0,T;H^{1}_{\sharp}(\Omega))}+\|\partial_{t}w_{n}\|_{L^{2}(0,T;H^{-1}_{\sharp}(\Omega))}\leq
C_{\vec{B},J}(\|w_{0}\|_{2}+\|v_{0}\|_{2}+\|F\|_{\infty}\sqrt{T}),$

and

$\|v_{n}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|v_{n}\|_{L^{2}(0,T;H^{1}_{\sharp}(\Omega))}+\|\partial_{t}v_{n}\|_{L^{2}(0,T;H^{-1}_{\sharp}(\Omega))}\leq
C_{\vec{B},J}(\|w_{0}\|_{2}+\|v_{0}\|_{2}+\|F\|_{\infty}\sqrt{T}),$

for each $n=1,2,\dots.$

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}\geq 0$, we may use the Lipschitz continuity of $F$.

######  Lemma 2.2.

There exists a $t^{*}>0$ such that the sequence $(w_{n},v_{n})$ is Cauchy in
$L^{\infty}(0,t^{*};L^{2}(\Omega))$.

###### Proof.

Let $\bar{w}_{n}:=w_{n}-w_{n-1}$ and $\bar{v}_{n}:=v_{n}-v_{n-1}$. Then from
(14) we have that $\bar{w}_{n+1}$ and $\bar{v}_{n+1}$ satisfy the weak form:

$\left\\{\begin{split}\langle\partial_{t}\bar{w}_{n+1},\psi\rangle+\int_{\Omega}\nabla\bar{w}_{n+1}\cdot\nabla\psi+w_{n}B_{1}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\cdot\nabla\psi+\bar{w}_{n+1}B_{1}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\psi\mathop{}\\!\mathrm{d}x\\\
=\int_{\Omega}[F(\tfrac{1}{2}(w_{n}+v_{n}))-F(\tfrac{1}{2}(w_{n-1}+v_{n-1}))]\psi\mathop{}\\!\mathrm{d}{x},\\\
\langle\partial_{t}\bar{v}_{n+1},\eta\rangle+\int_{\Omega}\nabla\bar{v}_{n+1}\cdot\nabla\eta+v_{n}B_{2}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\cdot\nabla\eta+\bar{v}_{n+1}B_{2}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\eta\mathop{}\\!\mathrm{d}x\\\
\quad=\int_{\Omega}[F(\tfrac{1}{2}(w_{n}+v_{n}))-F(\tfrac{1}{2}(w_{n-1}+v_{n-1}))]\eta\mathop{}\\!\mathrm{d}{x}.\end{split}\right.$
(15)

The above formulation allows us by testing with
$(\psi,\eta)=(\bar{w}_{n+1},\bar{v}_{n+1})$ to arrive at the identities

$\left\\{\begin{split}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\bar{w}_{n+1}\|_{2}^{2}+\|\nabla\bar{w}_{n+1}\|_{2}^{2}&=-\int_{\Omega}w_{n}B_{1}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\cdot\nabla\bar{w}_{n+1}+\bar{w}_{n+1}B_{1}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\bar{w}_{n+1}\mathop{}\\!\mathrm{d}x\\\
&\qquad+\int_{\Omega}[F(\tfrac{1}{2}(w_{n}+v_{n}))-F(\tfrac{1}{2}(w_{n-1}+v_{n-1}))]\bar{w}_{n+1}\mathop{}\\!\mathrm{d}{x},\\\
\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\bar{v}_{n+1}\|_{2}^{2}+\|\nabla\bar{v}_{n+1}\|_{2}^{2}&=-\int_{\Omega}v_{n}B_{2}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\cdot\nabla\bar{v}_{n+1}+\bar{v}_{n+1}B_{2}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\bar{v}_{n+1}\mathop{}\\!\mathrm{d}x\\\
&\qquad+\int_{\Omega}[F(\tfrac{1}{2}(w_{n}+v_{n}))-F(\tfrac{1}{2}(w_{n-1}+v_{n-1}))]\bar{v}_{n+1}\mathop{}\\!\mathrm{d}{x}.\end{split}\right.$
(16)

Working with the right-hand side of the first equation above we have the first
term,

$\displaystyle\int_{\Omega}w_{n}B_{1}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\cdot\nabla\bar{w}_{n+1}\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq\|B_{1}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\|_{\infty}\int_{\Omega}|w_{n}\nabla\bar{w}_{n+1}|\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq
C_{\vec{B}}C_{J}(\|\bar{w}_{n}\|_{2}+\|\bar{v}_{n}\|_{2})\,\|w_{n}\|_{2}\,\|\nabla\bar{w}_{n+1}\|_{2},$

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

$\displaystyle\int_{\Omega}w_{n}B_{1}(\nabla
J*(\bar{w}_{n}-\bar{v}_{n}))\cdot\nabla\bar{w}_{n+1}\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq
C_{\vec{B},J,F,w_{0},v_{0}}(\|\bar{w}_{n}\|_{2}^{2}+\|\bar{v}_{n}\|_{2}^{2})+\frac{1}{2}\|\nabla\bar{w}_{n+1}\|_{2}^{2}.$

Looking now to the second term, we have via Young’s inequality for
convolutions, Hölder’s inequality, and Young’s product inequality

$\displaystyle\int_{\Omega}\bar{w}_{n+1}B_{1}(\nabla
J*(w_{n}-v_{n}))\cdot\nabla\bar{w}_{n+1}\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq\|\nabla
J*(w_{n}-v_{n})\|_{\infty}\int_{\Omega}|\bar{w}_{n+1}\nabla\bar{w}_{n+1}|\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq
C_{J,F,w_{0},v_{0}}\|\bar{w}_{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla\bar{w}_{n+1}\|_{2}^{2}.$

Finally, from the third term on the right-hand side we see by Lipschitz
continuity of $F$ (denoting the Lipschitz constant by $L_{F}$),

$\displaystyle\int_{\Omega}[F(\tfrac{1}{2}(w_{n}+v_{n}))-F(\tfrac{1}{2}(w_{n-1}+v_{n-1}))]\bar{w}_{n+1}\mathop{}\\!\mathrm{d}{x}$
$\displaystyle\leq
L_{F}\int_{\Omega}\tfrac{1}{2}|\bar{w}_{n}+\bar{v}_{n}||\bar{w}_{n}|\mathop{}\\!\mathrm{d}{x}$
$\displaystyle\leq
C_{F}(\|\bar{w}_{n}\|_{2}^{2}+\|\bar{v}_{n}\|_{2}^{2}+\|\bar{w}_{n+1}\|_{2}^{2}),$

where in the last step we again make use of Hölder’s inequality and Young’s
inequality.

Repeating these calculations for $\bar{v}_{n+1}$ and summing the equations in
(16), we arrive at the estimates

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\bar{w}_{n+1}\|_{2}^{2}+\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\bar{v}_{n+1}\|_{2}^{2}\leq
C_{\vec{B},J,F,w_{0},v_{0}}(\|\bar{w}_{n}\|_{2}^{2}+\|\bar{v}_{n}\|_{2}^{2}+\|\bar{w}_{n+1}\|_{2}^{2}+\|\bar{v}_{n+1}\|_{2}^{2}).$
(17)

Choosing $t^{*}$ such that

$\frac{C_{\vec{B},J,F,w_{0},v_{0}}t^{*}}{1-C_{\vec{B},J,F,w_{0},v_{0}}t^{*}}<1$

and integrating (17) from 0 to $t<t^{*}$ we can arrive at the estimate

$\|\bar{w}_{n+1}(t,\cdot)\|_{2}^{2}+\|\bar{v}_{n+1}(t,\cdot)\|_{2}^{2}\leq
C_{\vec{B},J,F,w_{0},v_{0}}\left(\int_{0}^{t}\|\bar{w}_{n}(s,\cdot)\|_{2}^{2}+\|\bar{v}_{n}(s,\cdot)\|_{2}^{2}\mathop{}\\!\mathrm{d}s+\int_{0}^{t}\|\bar{w}_{n+1}(s,\cdot)\|_{2}^{2}+\|\bar{v}_{n+1}(s,\cdot)\|_{2}^{2}\mathop{}\\!\mathrm{d}s\right).$

Taking the supremum of both sides over $[0,t^{*})$, we obtain:

$\displaystyle\|\bar{w}_{n+1}\|^{2}_{L^{\infty}(0,t^{*};L^{2}(\Omega))}+\|\bar{v}_{n+1}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}$
$\displaystyle\leq
C_{\vec{B},J,F,w_{0},v_{0}}\left(\|\bar{w}_{n}\|_{L^{2}(0,t^{*};L^{2}(\Omega))}^{2}+\|\bar{v}_{n}\|_{L^{2}(0,t^{*};L^{2}(\Omega))}^{2}\right.$
$\displaystyle\quad+\left.\|\bar{w}_{n+1}\|_{L^{2}(0,t^{*};L^{2}(\Omega))}^{2}+\|\bar{v}_{n+1}\|_{L^{2}(0,t^{*};L^{2}(\Omega))}^{2}\right)$
$\displaystyle\leq
C_{\vec{B},J,F,w_{0},v_{0}}t^{*}\left(\|\bar{w}_{n}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}+\|\bar{v}_{n}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}\right.$
$\displaystyle\quad+\left.\|\bar{w}_{n+1}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}+\|\bar{v}_{n+1}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}\right),$

where we have used the inequality

$\|u\|_{L^{2}(0,t^{*};L^{2}(\Omega))}^{2}\leq
t^{*}\|u\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}.$ (18)

Manipulation of the above inequality yields

$\|\bar{w}_{n+1}\|^{2}_{L^{\infty}(0,t^{*};L^{2}(\Omega))}+\|\bar{v}_{n+1}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}\leq\frac{C_{\vec{B},J,F,w_{0},v_{0}}t^{*}}{1-C_{\vec{B},J,F,w_{0},v_{0}}t^{*}}\left(\|\bar{w}_{n}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}+\|\bar{v}_{n}\|_{L^{\infty}(0,t^{*};L^{2}(\Omega))}^{2}\right),$
(19)

which implies the sequence is Cauchy in $L^{\infty}(0,t^{*};L^{2}(\Omega))$.
Moreover, using inequality (18), we see this implies the sequence is also
Cauchy in $L^{2}(0,t^{*};L^{2}(\Omega))$. ∎

######  Remark.

The uniqueness of solutions to (12) follows from similar arguments presented
in the above proof. Indeed, letting $\bar{w}=w_{1}-w_{2}$ and
$\bar{v}=v_{1}-v_{2}$ where $(w_{1},v_{1})$ and $(w_{2},v_{2})$ are solutions
to (12), we see that $\bar{w}$ and $\bar{v}$ satisfy similar equations to
(15), namely

$\left\\{\begin{split}\langle\partial_{t}\bar{w},\psi\rangle+\int_{\Omega}\nabla\bar{w}\cdot\nabla\psi+w_{2}B_{1}(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla\psi+\bar{w}B_{1}(\nabla
J*(w_{1}-v_{1}))\cdot\nabla\psi\mathop{}\\!\mathrm{d}x\\\
=\int_{\Omega}[F(\tfrac{1}{2}(w_{1}+v_{1}))-F(\tfrac{1}{2}(w_{2}+v_{2}))]\psi\mathop{}\\!\mathrm{d}{x},\\\
\langle\partial_{t}\bar{v},\eta\rangle+\int_{\Omega}\nabla\bar{v}\cdot\nabla\eta+v_{2}B_{2}(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla\eta+\bar{v}B_{2}(\nabla
J*(w_{1}-v_{1}))\cdot\nabla\eta\mathop{}\\!\mathrm{d}x\\\
=\int_{\Omega}[F(\tfrac{1}{2}(w_{1}+v_{1}))-F(\tfrac{1}{2}(w_{2}+v_{2}))]\eta\mathop{}\\!\mathrm{d}{x}.\end{split}\right.$

Following the calculations which lead to equation (17), we arrive at the
differential inequality

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\bar{w}\|_{2}^{2}+\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|\bar{v}\|_{2}^{2}\leq
C_{\vec{B},J,F,w_{0},v_{0}}(\|\bar{w}\|_{2}^{2}+\|\bar{v}\|_{2}^{2}).$

A standard application of Grönwall’s inequality implies
$\|\bar{w}\|_{2}^{2}+\|\bar{v}\|_{2}^{2}=0$ for all $t\in[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^{\infty}(0,T;L^{2}(\Omega))^{2}\text{ (or
}L^{2}(Q_{T})^{2})$ and, from the energy estimates, is weakly convergent
(along a subsequence) in $L^{2}(0,T;H^{1}_{\sharp}(\Omega))$ with
$(\partial_{t}w_{n},\partial_{t}v_{n})$ weakly convergent (along a
subesquence) in $L^{2}(0,T;H^{-1}_{\sharp}(\Omega))$. Using these results, we
can pass to the limit in the weak formulation (14) (following the calculations
used to derive equation (17)) as $n\longrightarrow\infty$ to arrive at the
weak formulation of Definition 2.1. Thus, the limit is a solution to (12).

##  3 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:

$\left\\{\begin{split}\partial_{t}w&=\Delta w+\text{div}[\beta(m-1)(\nabla
J*(w-v))w]+F(\tfrac{1}{2}(w+v)),\\\ \partial_{t}v&=\Delta
v+\text{div}[\beta(m+1)(\nabla J*(w-v))v]+F(\tfrac{1}{2}(w+v)),\\\
m&=\frac{1}{2}(w-v),\end{split}\right.\quad(t,x)\in Q_{T}.$ (9)

Taking into account assumption (10), solutions to (9) are understood as in
Definition 2.1 with $B_{1}:=\beta(m-1)$ and $B_{2}:=\beta(m+1)$. We aim to
show the existence of solutions to a fixed-point to the map
$\mathcal{T}:\mathcal{Z}\longrightarrow\mathcal{Z}$ given by

$\bar{m}\xmapsto{\mathcal{T}_{1}}(\bar{w},\bar{v})\xmapsto{\mathcal{T}_{2}}\frac{1}{2}(\bar{w}-\bar{v}),$

where

$\mathcal{Z}:=\\{u\in L^{2}(Q_{T}):-1\leq u\leq 1\text{ almost everywhere in
}Q_{T}\\},$

$\mathcal{T}_{1}(\bar{m}):=(\bar{w},\bar{v})$ is the weak solution to (12)
with $B_{1}=\beta(\bar{m}-1)$ and $B_{2}=\beta(\bar{m}+1)$, and
$\mathcal{T}_{2}(\bar{w},\bar{v})=\frac{1}{2}(\bar{w}-\bar{v})$.

Now, using the work done in the previous sections, we have the following
lemmas:

######  Lemma 3.1.

The map $\mathcal{T}=\mathcal{T}_{2}\circ\mathcal{T}_{1}$ is well-defined over
the set $\mathcal{Z}$ and $\mathcal{T}:\mathcal{Z}\longrightarrow\mathcal{Z}$.

###### Proof.

We study separately the maps $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$.

#### The map $\mathcal{T}_{1}$ :

Fix $\bar{m}\in\mathcal{Z}$. By Theorem 2.1, we have the existence of a unique
pair

$(\bar{w},\bar{v})\in\mathcal{X}^{2}\text{ with
}(\partial_{t}\bar{w},\partial_{t}\bar{v})\in\left(L^{2}(0,T;H^{-1}_{\sharp}(\Omega))\right)^{2}$

such that $0\leq\bar{w},\bar{v}$. Thus the map $\mathcal{T}_{1}$ is well-
defined.

#### The map $\mathcal{T}_{2}$ :

The map $\mathcal{T}_{2}$ is trivially well-defined. What remains to be shown
is that $\mathcal{T}_{2}$ maps into $\mathcal{Z}$. To see this, define
$m:=\frac{1}{2}(\bar{w}-\bar{v})\in L^{2}(Q_{T})$ and notice that

$\partial_{t}m=\Delta
m-\text{div}[\beta(\bar{w}+\bar{v}-\bar{m}(\bar{w}-\bar{v}))\nabla J*m],$ (20)

with $|m_{0}|\leq 1.$ Since $\bar{m}\in L^{\infty}(Q_{T})$ and
$\bar{w},\bar{v}\in L^{2}(0,T;H^{1}_{\sharp}(\Omega))\cap
L^{\infty}(0,T;L^{2}(\Omega))$, we have
$\beta(\bar{w}+\bar{v}-\bar{m}(\bar{w}-\bar{v}))=:B(t,x)\in
L^{\infty}(0,T;L^{2}(\Omega))$ . Thus, by Proposition 2.2, we have that
$|m|\leq 1$. ∎

######  Remark 3.1.

With similar arguments to those found in Proposition 2.1, it can be easily
shown the function $\phi:=\frac{1}{2}(\bar{w}+\bar{v})$, which satisfies the
equation

$\partial_{t}\phi=\Delta\phi-\text{div}[2\beta\bar{m}(1-\phi)\nabla
J*m]+F(\phi),$

is nonnegative and bounded by 1 (i.e.,$0\leq\phi\leq 1$) since
$\bar{m},m\in\mathcal{Z}$ and $0\leq\bar{w},\bar{v}$. Together, these bounds
allow us to conclude $0\leq\bar{w},\bar{v}\leq 2$ almost everywhere in
$Q_{T}$.

######  Lemma 3.2.

The map
$\mathcal{T}=\mathcal{T}_{2}\circ\mathcal{T}_{1}:\mathcal{Z}\longrightarrow\mathcal{Z}$
is continuous.

###### Proof.

Let $\bar{m}_{n},\bar{m}\in\mathcal{Z}$ be such that
$\bar{m}_{n}\longrightarrow\bar{m}$ in $L^{2}(Q_{T})$. Define also
$(\bar{w}_{n},\bar{v}_{n}):=\mathcal{T}_{1}(m_{n})$ and
$(\bar{w},\bar{v}):=\mathcal{T}_{1}(m).$ We want to show
$(\bar{w}_{n},\bar{v}_{n})\longrightarrow(\bar{w},\bar{v})$ in $L^{2}(Q_{T}).$
To this end, let $m:=\bar{m}_{n}-\bar{m}$, $w:=\bar{w}_{n}-\bar{w}$, and
$v:=\bar{v}_{n}-\bar{v}$ for fixed $n$. Following similar manipulations which
gave rise to equation (15), we have that $w$ and $v$ satisfy

$\left\\{\begin{split}\langle\partial_{t}w,\psi\rangle+\int_{\Omega}\nabla
w\cdot\nabla\psi\mathop{}\\!\mathrm{d}x&=-\int_{\Omega}\beta\bar{w}_{n}(\bar{m}_{n}-1)(\nabla
J*(w-v))\cdot\nabla\psi+w(\bar{m}_{n}-1)(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla\psi\mathop{}\\!\mathrm{d}x\\\
&\qquad-\int_{\Omega}\bar{w}\,m(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla\psi\mathop{}\\!\mathrm{d}x+\int_{\Omega}[F(\tfrac{1}{2}(\bar{w}_{n}+\bar{v}_{n})-F(\tfrac{1}{2}(\bar{w}+\bar{v}))]\psi\mathop{}\\!\mathrm{d}{x},\\\
\langle\partial_{t}v,\eta\rangle+\int_{\Omega}\nabla
v\cdot\nabla\eta\mathop{}\\!\mathrm{d}x&=-\int_{\Omega}\beta\bar{v}_{n}(\bar{m}_{n}+1)(\nabla
J*(w-v))\cdot\nabla\eta+v(\bar{m}_{n}+1)(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla\eta\mathop{}\\!\mathrm{d}x\\\
&\qquad-\int_{\Omega}\bar{v}\,m(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla\eta\mathop{}\\!\mathrm{d}x+\int_{\Omega}[F(\tfrac{1}{2}(\bar{w}_{n}+\bar{v}_{n})-F(\tfrac{1}{2}(\bar{w}+\bar{v}))]\eta\mathop{}\\!\mathrm{d}{x}.\end{split}\right.$
(21)

Testing with $(\psi,\eta)=(w,v)$ we have the following equations:

$\left\\{\begin{split}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|w\|_{2}^{2}+\|\nabla
w\|_{2}^{2}&=-\int_{\Omega}\beta\bar{w}_{n}(\bar{m}_{n}-1)(\nabla
J*(w-v))\cdot\nabla w+w(\bar{m}_{n}-1)(\nabla J*(\bar{w}-\bar{v}))\cdot\nabla
w\mathop{}\\!\mathrm{d}x\\\ &\qquad-\int_{\Omega}\bar{w}\,m(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla
w\mathop{}\\!\mathrm{d}x+\int_{\Omega}[F(\tfrac{1}{2}(\bar{w}_{n}+\bar{v}_{n})-F(\tfrac{1}{2}(\bar{w}+\bar{v}))]w\mathop{}\\!\mathrm{d}{x},\\\
\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|v\|_{2}^{2}+\|\nabla
v\|_{2}^{2}&=-\int_{\Omega}\beta\bar{v}_{n}(\bar{m}_{n}+1)(\nabla
J*(w-v))\cdot\nabla v+v(\bar{m}_{n}+1)(\nabla J*(\bar{w}-\bar{v}))\cdot\nabla
v\mathop{}\\!\mathrm{d}x\\\ &\qquad-\int_{\Omega}\bar{v}\,m(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla
v\mathop{}\\!\mathrm{d}x+\int_{\Omega}[F(\tfrac{1}{2}(\bar{w}_{n}+\bar{v}_{n})-F(\tfrac{1}{2}(\bar{w}+\bar{v}))]v\mathop{}\\!\mathrm{d}{x}.\end{split}\right.$
(22)

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:

$\left\\{\begin{split}\left|\int_{\Omega}\beta\left[\bar{w}_{n}(\bar{m}_{n}-1)(\nabla
J*(w-v))+w(\bar{m}_{n}-1)(\nabla J*(\bar{w}-\bar{v}))\right]\cdot\nabla
w\mathop{}\\!\mathrm{d}x\right|&\leq
C_{J,F}(\|w\|_{2}^{2}+\|v\|_{2}^{2})+\frac{1}{3}\|\nabla w\|_{2}^{2},\\\
\left|\int_{\Omega}\beta\left[\bar{v}_{n}(\bar{m}_{n}+1)(\nabla
J*(w-v))+v(\bar{m}_{n}+1)(\nabla J*(\bar{w}-\bar{v}))\right]\cdot\nabla
v\mathop{}\\!\mathrm{d}x\right|&\leq
C_{J,F}(\|w\|_{2}^{2}+\|v\|_{2}^{2})+\frac{1}{3}\|\nabla
v\|_{2}^{2}.\end{split}\right.$ (23)

Similarly, the evaporation terms follow the same arguments,

$\left\\{\begin{split}\left|\int_{\Omega}[F(\tfrac{1}{2}(\bar{w}_{n}+\bar{v}_{n})-F(\tfrac{1}{2}(\bar{w}+\bar{v}))]w\mathop{}\\!\mathrm{d}x\right|&\leq
C_{F}(\|w\|_{2}^{2}+\|v\|_{2}^{2}),\\\
\left|\int_{\Omega}[F(\tfrac{1}{2}(\bar{w}_{n}+\bar{v}_{n})-F(\tfrac{1}{2}(\bar{w}+\bar{v}))]v\mathop{}\\!\mathrm{d}x\right|&\leq
C_{F}(\|w\|_{2}^{2}+\|v\|_{2}^{2}).\end{split}\right.$ (24)

Meanwhile, the new third term yields the bound

$\left\\{\begin{split}\left|\int_{\Omega}\bar{w}\,m(\nabla
J*(\bar{w}-\bar{v}))\cdot\nabla w\mathop{}\\!\mathrm{d}x\right|&\leq
C_{J,F}\|m\|_{2}^{2}+\frac{1}{3}\|\nabla w\|_{2}^{2},\\\
\left|\int_{\Omega}\beta\bar{v}\,m(\nabla J*(\bar{w}-\bar{v}))\cdot\nabla
v\mathop{}\\!\mathrm{d}x\right|&\leq C_{J,F}\|m\|_{2}^{2}+\frac{1}{3}\|\nabla
v\|_{2}^{2}.\end{split}\right.$ (25)

Summing up the two inequalities and maximizing constants leaves us with

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|w\|_{2}^{2}+\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|v\|_{2}^{2}\leq
C_{J,F}(\|m\|_{2}^{2}+\|w\|_{2}^{2}+\|v\|_{2}^{2}).$ (26)

Integrating the above equation from $0$ to $t$ gives us the integral
inequality

$\|w\|_{2}^{2}+\|v\|_{2}^{2}\leq
C_{J,F}\|m\|_{L^{2}((0,t)\times\Omega)}^{2}+C_{J,F}\int_{0}^{t}\|w\|_{2}^{2}+\|v\|_{2}^{2}\mathop{}\\!\mathrm{d}s.$
(27)

Grönwall’s inequality and integration provides us with the estimate:

$\|w\|_{L^{2}(Q_{T})}^{2}+\|v\|_{L^{2}(Q_{T})}^{2}\leq
C_{J,F}T\|m\|_{L^{2}(Q_{T})}^{2}\exp{(C_{J,F}T)}.$ (28)

Continuity of $\mathcal{T}_{1}$ and therefore $\mathcal{T}$ follows. ∎

######  Theorem 3.1.

The map $\mathcal{T}:=\mathcal{T}_{2}\circ\mathcal{T}_{1}$ has at least one
fixed-point $m^{*}\in\mathcal{W}\cap\mathcal{Z}.$

###### Proof.

Summing up the results from before, we have the set $\mathcal{Z}$ is a non-
empty, closed, and convex subset of $L^{2}(Q_{T})$, the map
$\mathcal{T}:\mathcal{Z}\longrightarrow\mathcal{Z}$ is well-defined and
continuous (Lemmas 3.1 and 3.2), and
$\mathcal{T}(\mathcal{Z})\subseteq\mathcal{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
$\mathcal{W}\cap\mathcal{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=\frac{1}{2}(w-v)$ and $\phi=\frac{1}{2}(w+v)$ are solutions to equation
(7). Moreover, we have shown that $0\leq w,v\leq 2$ which implies that

$|m|\leq\phi.$

Since in the fixed-point argument we show that $m\in\mathcal{Z}$, we can
conclude via Proposition 2.1 that $\phi\leq 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=\tfrac{1}{2}(w-v)$, we have that inequality (26) now reads

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|w\|_{2}^{2}+\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\|v\|_{2}^{2}\leq
C_{J,F}(\|w\|_{2}^{2}+\|v\|_{2}^{2}).$

Grönwall’s inequality then proves the result. ∎

This proves Theorem 1.1.

##  4 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.

###  4.1 Simulation results

To fix ideas, we take into account a linear evaporation model, i.e.,
$F:\mathbb{R}\to\mathbb{R}$ defined by $F(r)=\alpha(1-r)$ for all
$r\in\mathbb{R}$ with given $\alpha>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 $\Delta t,\Delta x,\Delta y>0$ we discretize the domain $\Omega$ by
cells
$\Lambda_{i,j}:=[x_{i}-\frac{1}{2}h,x_{i}+\frac{1}{2}h)\times[y_{j}-\frac{1}{2}h,y_{j}+\frac{1}{2}h)$
for $i,j=1,2,\dots,N$. To account for the periodic boundary conditions, we
periodically extend functions defined on the nodes, i.e., $f_{i,j}=f_{i\pm
N,j}=f_{i,j\pm N}$. The initial pair $(m_{0},\phi_{0})$ is then approximated
by

$m^{0}_{i,j}:=\frac{1}{|\Lambda_{i,j}|}\int_{\Lambda_{i,j}}m_{0}(x,y)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y\quad\text{
and
}\quad\phi^{0}_{i,j}:=\frac{1}{|\Lambda_{i,j}|}\int_{\Lambda_{i,j}}\phi_{0}(x,y)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$

and define the fully explicit scheme:

$\left\\{\begin{split}\frac{m^{k+1}_{i,j}-m^{k}_{i,j}}{\Delta
t}&=\frac{1}{\Delta x^{2}}D^{2}_{i}[m^{k}_{i,j}]+\frac{1}{\Delta
y^{2}}D^{2}_{j}[m^{k}_{i,j}]\\\ &-\frac{\beta}{\Delta
x}D^{1}_{i}[(\phi^{k}_{i,j}-(m^{k}_{i,j})^{2})\tilde{J}^{k}_{x,i,j}]-\frac{\beta}{\Delta
y}D^{1}_{j}[(\phi^{k}_{i,j}-(m^{k}_{i,j})^{2})\tilde{J}^{k}_{y,i,j}]\\\
\frac{\phi^{k+1}_{i,j}-\phi^{k}_{i,j}}{\Delta t}&=\frac{1}{\Delta
x^{2}}D^{2}_{i}[\phi^{k}_{i,j}]+\frac{1}{\Delta
y^{2}}D^{2}_{j}[\phi^{k}_{i,j}]\\\ &-\frac{\beta}{\Delta
x}D^{1}_{i}[m^{k}_{i,j}(1-\phi^{k}_{i,j})\tilde{J}^{k}_{x,i,j}]-\frac{\beta}{\Delta
y}D^{1}_{j}[m^{k}_{i,j}(1-\phi^{k}_{i,j})\tilde{J}^{k}_{y,i,j}]+F(\phi_{i,j}^{k})\end{split}\right.,$

where

$D^{2}_{l}[f_{l}]:=f_{l+1}-2f_{l}+f_{l-1},\quad
D^{1}_{l}[f_{l}]:=f_{l+1}-f_{l-1},$

and $\tilde{J}^{k}_{x,i,j}$ (respectively, $\tilde{J}^{k}_{y,i,j}$) denotes
the approximation of $\partial_{x}J*m^{k}(x_{i},y_{j})$ (respectively,
$\partial_{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(\phi)\equiv 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.

Figure 1: Simulation of equation (7) with $\alpha=0.1$, $\beta=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 $\phi$ is near one are colored blue whereas
regions near zero are colored red.

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-\phi\|_{1}$. As expected, we observe a decrease in solvent ratio over
time. Next, we plot the $\|\cdot\|_{1}$ norms of both $m$ and $\phi$ (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.

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.

##  5 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|-\phi|\longrightarrow 0$ almost everywhere as
$t\longrightarrow\infty$ 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.

## Acknowledgments

The authors are grateful to the anonymous reviewers for their comments which
significantly contributed to the quality of the manuscript. The authors also
thank M. Eden (Karlstad, Sweden) for fruitful discussions on the mathematical
analysis of the auxiliary problem and S. A. Muntean (Karlstad) for
brainstorming about the modeling of the evaporation production term. RL and AM
acknowledge the financial support of Carl Tryggers Stieftelse via the grant
CTS 21:1656.

## References

  * [1] H. Amann.  Ordinary differential equations: an introduction to nonlinear analysis, volume 13.  Walter de gruyter, 1990.
  * [2] A. K. Barua, R. Chew, S. Li, J. Lowengrub, A. Münch, and B. Wagner.  Sharp-interface problem of the Ohta-Kawasaki model for symmetric diblock copolymers.  Journal of Computational Physics, 2023.
  * [3] J. A. Carrillo, Y.-P. Choi, and M. Hauray.  The derivation of swarming models: mean-field limit and Wasserstein distances.  In A. Muntean and F. Toschi, editors, Collective Dynamics from Bacteria to Crowds, pages 1–46. CISM Series, 2014.
  * [4] J. A. Carrillo, R. S. Gvalani, G. A. Pavliotis, and A. Schlichting.  Long-time behaviour and phase transitions for the McKean–Vlasov equation on the torus.  Archive for Rational Mechanics and Analysis, 235(1):635–690, 2020\.
  * [5] B. Chazelle, Q. Jiu, Q. Li, and C. Wang.  Well-posedness of the limiting equation of a noisy consensus model in opinion dynamics.  Journal of Differential Equations, 263(1):365–397, 2017.
  * [6] E. N. M. Cirillo, M. Colangeli, E. Moons, A. Muntean, S. A. Muntean, and J. van Stam.  A lattice model approach to the morphology formation from ternary mixtures during the evaporation of one component.  Eur. Phys. J. Spec. Top., 228:55–68, 2019.
  * [7] C. Creton and M. Ciccotti.  Fracture and adhesion of soft materials.  Reports on Progress in Physics, 79(4):046601, 2016.
  * [8] J. Cummings, J. S. Lowengrub, B. G. Sumpter, S. M. Wise, and R. Kumar.  Modeling solvent evaporation during thin film formation in phase separating polymer mixtures.  Soft Matter, 14(10):1833–1846, 2018.
  * [9] M. Eden, C. Nikolopoulos, and A. Muntean.  A multiscale quasilinear system for colloids deposition in porous media: Weak solvability and numerical simulation of a near-clogging scenario.  Nonlinear Analysis: Real World Applications, 63:103408, 2022.
  * [10] L. C. Evans.  Partial Differential Equations, volume 19.  American Mathematical Society, $2^{nd}$ edition, 2010.
  * [11] C. G. Gal, A. Giorgini, and M. Grasselli.  The nonlocal Cahn–Hilliard equation with singular potential: well-posedness, regularity and strict separation property.  Journal of Differential Equations, 263(9):5253–5297, 2017.
  * [12] C. G. Gal, A. Giorgini, and M. Grasselli.  The separation property for 2D Cahn-Hilliard equations: local, nonlocal and fractional energy cases.  Discrete and Continuous Dynamical Systems, 10, 2022.
  * [13] C. G. Gal and M. Grasselli.  Longtime behavior of nonlocal Cahn-Hilliard equations.  Discrete & Continuous Dynamical Systems, 34(1):145, 2014.
  * [14] G. Giacomin and J. L. Lebowitz.  Phase segregation dynamics in particle systems with long range interaction I. Macroscopic limits.  J. Statist. Phys., 87:37–61, 1997.
  * [15] G. Giacomin and J. L. Lebowitz.  Phase segregation dynamics in particle systems with long range interactions II: Interface motion.  SIAM Journal on Applied Mathematics, 58(6):1707–1729, 1998.
  * [16] Z. Guan, J. S. Lowengrub, C. Wang, and S. M. Wise.  Second order convex splitting schemes for periodic nonlocal Cahn–Hilliard and Allen–Cahn equations.  Journal of Computational Physics, 277:48–71, Nov 2014.
  * [17] Z. Guan, C. Wang, and S. M. Wise.  A convergent convex splitting scheme for the periodic nonlocal Cahn-Hilliard equation.  Numerische Mathematik, 128(2):377–406, Oct 2014.
  * [18] M. G. Hennessy, V. M. Burlakov, A. Goriely, B. Wagner, and A. Münch.  Controlled topological transitions in thin-film phase separation.  SIAM Journal on Applied Mathematics, 75(1):38–60, 2015.
  * [19] H. Hoppe and N. S. Sariciftci.  Organic solar cells: An overview.  Journal of Materials Research, 19(7):1924–1945, 2004.
  * [20] V. A. Khoa and A. Muntean.  Corrector homogenization estimates for a non-stationary Stokes-Nernst-Planck-Poisson system in perforated domains.  Communications in Mathematical Sciences, 17(3):705–738, 2019.
  * [21] J.-M. Lasry and P.-L. Lions.  Mean field games.  Japanese Journal of Mathematics, 2(1):229–260, 2007.
  * [22] W. Liu, L. Wu, and C. Zhang.  Long-time behaviors of mean-field interacting particle systems related to McKean–Vlasov equations.  Communications in Mathematical Physics, 387(1):179–214, 2021.
  * [23] R. Lyons, S. A. Muntean, E. N. M. Cirillo, and A. Muntean.  A continuum model for morphology formation from interacting ternary mixtures: Simulation study of the formation and growth of patterns.  Physica D: Nonlinear Phenomena, 453:133832, 2023.
  * [24] R. Marra and M. Mourragui.  Phase segregation dynamics for the Blume–Capel model with Kac interaction.  Stochastic Processes and their Applications, 88(1):79–124, 2000\.
  * [25] A. Miranville.  The Cahn-Hilliard Equation: Recent Advances and Applications.  CBMS-NSF Regional Conference Series in Applied Mathematics. SIAM, 2019\.
  * [26] M. Müller, A. Lang, M. Klüppel, and U. Giese.  Influence of phase morphology on viscoelastic properties of rubber blends.  In C. Marano, F. B. Vangosa, L. Andena, and R. Frassine, editors, Constitutive Models for Rubber XII: Proceedings of the 12th European Conference on Constitutive Models for Rubber (ECCMR 2022). CRC Press, 2022.
  * [27] A. Muntean.  Continuum Modeling: An Approach Through Practical Examples.  SpringerBriefs in Mathematical Methods. Springer International Publishing, 2015.
  * [28] S. A. Muntean, V. C. E. Kronberg, M. Colangeli, A. Muntean, J. van Stam, E. Moons, and E. N. M. Cirillo.  Quantitative analysis of phase formation and growth in ternary mixtures upon evaporation of one component.  Phys. Rev. E, 106:025306, 2022.
  * [29] J. Nelson, J. J. Kwiatkowski, J. Kirkpatrick, and J. M. Frost.  Modeling charge transport in organic photovoltaic materials.  Acc. Chem. Res., 42(11):1768–1778, 2009.
  * [30] R. L. Pego.  Front migration in the nonlinear Cahn-Hilliard equation.  Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 422(1863):261–278, 1989.
  * [31] E. Presutti.  Scaling Limits in Statistical Mechanics and Microstructures in Continuum Mechanics.  Theoretical and Mathematical Physics. Springer, 2008.
  * [32] R. Rabani, H. Sadafi, H. Machrafi, M. Abbasi, B. Haut, and P. Dauby.  Influence of evaporation on the morphology of a thin film of a partially miscible binary mixture.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 612:126001, 2021.
  * [33] C. Schaefer.  Theory of nanostructuring in solvent-deposited thin polymer films.  PhD thesis, Technische Universiteit Eindhoven, 2016.
  * [34] M. Setta, V. C. E. Kronberg, S. A. Muntean, E. Moons, J. van Stam, E. N. M. Cirillo, M. Colangeli, and A. Muntean.  A mesoscopic lattice model for morphology formation in ternary mixtures with evaporation.  Communications in Nonlinear Science and Numerical Simulation, 119:107083, 2023.
  * [35] A. K. Tiwari, A. Pandey, J. Paul, and A. Anand.  Fast accurate approximation of convolutions with weakly singular kernel and its applications, 2021.  (arXiv.2107.03958).
  * [36] J. M. R. Vera.  A convergent iterative method for a logistic chemotactic system.  Revista Colombiana de Matemáticas, 51(1):103–117, 2017.