A Mesoscopic Lattice Model for Morphology Formation in Ternary Mixtures with Evaporation

We develop a mesoscopic lattice model to study the morphology formation in interacting ternary mixtures with evaporation of one component. As concrete application of our model, we wish to capture morphologies as they are typically arising during fabrication of organic solar cells. In this context, we consider an evaporating solvent into which two other components are dissolved, as a model for a 2-component coating solution that is drying on a substrate. We propose a 3-spins dynamics to describe the evolution of the three interacting species. As main tool, we use a Monte Carlo Metropolis-based algorithm, with the possibility of varying the system's temperature, mixture composition, interaction strengths, and evaporation kinetics. The main novelty is the structure of the mesoscopic model -- a bi-dimensional lattice with periodic boundary conditions and divided in square cells to encode a mesoscopic range interaction among the units. We investigate the effect of the model parameters on the structure of the resulting morphologies. Finally, we compare the results obtained with the mesoscopic model with corresponding ones based on an analogous lattice model with a short range interaction among the units, i.e. when the mesoscopic length scale coincides with the microscopic length scale of the lattice.


Introduction
The mechanisms underlying mesoscopic and macroscopic pattern formation from local microscopic interactions are explored in many fields of physics, chemistry, and biology [18]. Lattice-based modelling of interactions between units 1 (magnetic spins, agents, molecules, pedestrians, colloids, etc.) can give a coherent description of real behavior in many different situations. Well-known examples include descriptions of phase transitions, flame propagation, spinodal decomposition, formation of magnetization bands, acceleration shock waves in traffic flow, building of coherent groups in large pedestrian crowds and molecules moving inside a cell searching for exit gates; see e.g. [13,24] and references cited therein. In this paper, we focus on a setting inspired from work done on organic solar cells (compare e.g. [12]), where the effect of the evaporation of a solvent -background environment for a mixture of two interacting polymers -on the formation of polymer-polymer stable mesoscopic configurations, called here morphologies, is of strong interest. From the energy harvesting perspective, this is a particularly relevant subject, since one expects that the shape and spatial arrangement of morphologies can affect considerably the overall power conversion efficiency of organic solar cells [6,19,33]. Conceptually related situations arise in the dynamics of interacting populations driven by different opinions and targets [23]. According to [24], both these applications belong to the realm of "complexity".
Using the framework offered by lattice-based modeling (see e.g. [1,20] for the general methodology and [10] for a review of the theoretical foundations), we aim at understanding to which extent the formation of stable spatial patterns (morphologies), which are obtained as a result of pair-wise interactions in a ternary mixture with one evaporating component, depend on specific length scales higher than the grid size of the lattice. We refer to such larger scales as mesoscales and we label them by λ. The interest in unveiling mesoscale effects was triggered by our previous simulation results reported in [4,15] and [21], where we noticed the occurrence of different types of morphology shapes. The simple observation that the geometry of the shapes depends on the choice of model parameters (the system's temperature, volatility, interaction parameters, etc.) makes us wonder whether our simulation results are going to be drastically different if the overall dynamics combines information not only from a microscopic scale, but from both microscopic (λ = 1) and mesoscopic length scales (1 < λ L, with L denoting some macroscopic length scale, e.g. the size of the simulation box). In particular, the vast separation of micro-, meso-and macro-scales is often invoked in the mathematical derivations of macroscopic behavior, used in statistical mechanics and in kinetic theory of gases as described for instance in [2,5,25].
Note that in [4] we explored the parameter space leading to morphologies as obtained with simulations done exclusively at the microscopic scale (i.e., pair-wise interactions involving only the nearest neighbours of randomly selected lattice sites). Further motivation towards performing mesoscale-level simulations in a closely related context is mentioned in [34].
In this paper, we propose a two-scale lattice model capable of producing morphologies, where the interactions within the mixture capture not only microscopic information (from pairs of spins) but also mesoscopic information (from pairs of λ-sized blocks of spins). The structure of such a two-scale model is inspired by the setting considered in Ch. 4 of [25], where one considers systems characterized by having the inter-atomic and mesoscopic characteristic interaction lengths sharply separated. In this context, the discussion is done in terms of Kac potentials; compare loc. cit.. Potentially, such an approach can provide an alternative two-scale model for which the so-called Lebowitz-Penrose mean-field limit might be proven rigorously, at least in the absence of the evaporation process, which poses additional mathematical challenges as it is a non-equilibrium interface process; see [17] for details on the Lebowitz-Penrose scaling of Hamiltonians and [25] for suitable mathematical techniques to study rigorously the passage to the continuum limit. We will study elsewhere the passage to the hydrodynamic limit in our setting. As we will see in Section 2, the proposed model incorporates both short-distance and long-distance interactions. Interestingly, a non-intuitive effect stands out -if a relatively low amount of solvent (evaporating component) is present in the mixture, then our two-scale model exhibits morphologies with self-similar features, i.e. one is able to zoom in and zoom out inside the geometry of the morphologies by suitably varying the intermediate scale λ together with a proportional modification of the size of the simulation box. This effect is a direct consequence of the scaling choice of the Hamiltonian functional driving the dynamics. On the other hand, with the current scaling we are unable to obtain mean-field effects. More work is needed in this direction, especially if one is tempted to find a rigorous link between the morphologies obtained with microscopic and/or mesoscopic lattice models, as done here, with the ones obtained by means of coupled systems of Cahn-Hilliard-type models, as for instance was performed in [26,30,32] and references cited therein.
The paper is organized as follows: In Section 2, we describe the proposed mesoscopic lattice model. The main work consists in performing simulation tests and interpreting the results. In Section 3, we show basic simulation results obtained for a fixed value of λ, bringing the attention to a particular mesoscopic level. To fix ideas, we have selected λ = 4. We are using this scenario to explore the effect of various parameters (mixture composition, temperature of the system, volatility, interaction strengths among mixture components) on the formation of morphologies. This type of numerical results show that our mesoscopic model is able to capture the type of results obtained in [4]. Additionally, a few typical mesoscopic features (like dependence of the morphology widths on λ and quicker stabilization of morphologies) can now be pointed out -such features arise at each allowed choice of λ > 1 and set the foundation for what we refer to as near self-similarity in morphology arrangements. The role of Section 4 is to show specific effects that can now be probed varying the mesoscale length λ. Here we discover ways to investigate the interplay between changes in the model parameters and size effects, as well as multiscale effects by observing the basic simulation output (check, e.g., Section 3) for distinct levels of λ. In the context of this section, we also study the connection between changing λ and varying the box size of the simulations. Section 5 concludes the paper with a discussion of our main findings as well as with an outlook on possible further research concerning this type of interacting mixtures and related matters.

Model description
To build our mesoscopic lattice model, we consider a two-dimensional rectangular lattice Λ := {1, . . . , L 1 }× {1, . . . , L 2 } with L 1 , L 2 ∈ N. We endow the lattice dynamics with periodic boundary conditions. An element of the torus Λ is called site. To reach mesoscopic descriptions, we introduce the long-distance interaction parameter λ. Concretely, we consider a partitioning of the lattice Λ in square cells of side λ min{L 1 , L 2 }. We end up with a lattice composed of l 1 × l 2 squares containing λ 2 sites each, where l 1 := L 1 /λ and l 2 := L 2 /λ. For computational convenience, we choose the values of λ, L 1 , and L 2 such that l 1 , l 2 ∈ N. Any cell X will be identified with the pair of integers (X 1 , X 2 ) ∈ {1, . . . , l 1 } × {1, . . . , l 2 }. Two cells X, Y such that X = Y are said to be nearest neighbors if their euclidean distance is one. We refer to the set of the nearest neighbours of the cell X as N (X). A pair of neighboring cells will be called bond.
Each site x ∈ Λ has an associated spin variable σ x ∈ {−1, 0, +1}. To model interaction we introduce the symmetric interaction tensor J ∈ R 3×3 and we adopt the notation J α,β with α, β = −1, 0, +1. We denote by σ ∈ {−1, 0, +1} Λ any configuration of the system on the lattice Λ and by σ ∆ its restriction to ∆ ⊂ Λ. In Figure 1, we show the subset ∆ = X ∪ Y , composed by the cells X, Y and their nearest neighbours, in different orientations. As main modelling step, given two neighboring cells X and Y and given two sites x ∈ X and y ∈ Y , we consider the local energy given by the Hamiltonian where for any α and Z the notation n Z α refers to the number of particles of the species α in the cell Z, see Figure 1.
In this case, the first two terms of the Hamiltonian capture the intra-cell energy, i.e. a mean field-like Hamiltonian restricted to the two cells X and Y (as in [25]) with respect to the fixed sites x and y. This is given by the sum of all the interactions between the spins at x and y with all the other sites in the corresponding cells, according to the interaction tensor J. Hence, the first two terms in (1) are the part of the Hamiltonian counting contributions from the microscopic scale. We involve a dimensionless tuning parameter C to control the weight of the intra-cell energy with respect to the total energy H x,y . The last two terms of the local Hamiltonian are the cell-cell interfacial energy, that represents the inter-cell energy. This part of the energy takes into account the interaction parameters of each combination of species weighted by the number of particles of those species in nearest neighbours cells and builds the mesoscopic contribution in the Hamiltonian. Since the Hamiltonian includes interactions at both microscopic and mesoscopic length scales, the resulting model is typically a multiscale model, or better said, a micro-meso model.
Usually, the interaction between particles of the same species is the one with minimum energy. Hence, to achieve the formation of morphologies, we want the energy computed with this Hamiltonian to be decreasing during the evolution of the system. We fix p −1 , p 0 , p +1 ∈ [0, 1) as initial probabilities to occupy one site such that p −1 + p 0 + p +1 = 1. We choose the initial configuration 2 σ 0 by setting any spin σ 0 x = i, i ∈ {−1, 0, +1} for all x ∈ Λ according to the probabilities p −1 , p 0 and p +1 . Sometimes, it will make sense to replace the probabilities p i with the actual percentages. We can now define p * := p +1 /(1 − p 0 ). We make use also of two additional parameters: the volatility 3 parameter φ ∈ [0, 1] and a parameter β > 0, such that k B T = β −1 is the thermal energy of the considered system, involving the Boltzmann constant k B . Finally, we fix as stop parameter, called p 0 , a certain percentage of the remaining solvent in the lattice. This stop parameter, called stop in the code, defines implicitly the final observation time T f of our mesoscopic model. It is worth noting that in our approach we cannot link directly the stop parameter with the correct estimation of the real time needed to build the morphology. Instead, to quantify the length of the simulations, we record the number of iterations needed to satisfy the stopping criterion. Therefore, the last configuration is denoted by σ T f and will be the one with a percentage of remaining solvent less than or equal to the one defined by the stop parameter p 0 . We consider the lattice to be oriented as in Figure 2. We divide the lattice in square cells of size λ min{L 1 , L 2 } to include the effects of a long distance interaction energy. The mathematical model is based on the discrete-time Markov chain Monte Carlo method [1]. We propose a Metropolis-type algorithm, completed with volatility and evaporation rules, as follows: 1. select a cell bond (X, Y ) in a uniformly random fashion; 2. select a site for each cell in the bond uniformly at random, namely x ∈ X and y ∈ Y ; 3. if X = (l 1 , X 2 ), Y = (1, X 2 ) and σ y = 0, then replace the zero with +1 with probability p * , configuration is called σ t . For convenience, the configuration at t = 0, i.e. σ 0 , is called "A" in the code, while for t > 0, σ t is called S. Additionally, we define a configuration S1 for the code. The latter is the configuration S with two swapped spins, as explained later on in this section. 3 The term "volatility" refers here to a drift leading the solvent particles towards the evaporation surface. It has nothing to do with the physical concept of volatility. The analogue name is though used since such break of symmetry applies only to one of the components of the mixture, i.e. to the potentially "volatile" one. otherwise with −1. Update p * according to the new percentages of particles p −1 , p 0 , p +1 . This is the evaporation of the zero (red) particle; 4. if X = (X 1 , X 2 ), Y = (X 1 + 1, X 2 ) and y = 0, with X 1 < l 1 , then exchange the two spins with probability φ. This is the volatility of the zero (red) particle; 5. otherwise swap the two spins, then consider the unswapped configuration S and the swapped one S1 on X ∪ Y . Compute ∆H = H x,y (S1) − H x,y (S) according to the Hamiltonian (1), that is the difference between the evolution of the Hamiltonian for the new and for the previous configuration. Accept the exchange with probability 1 if ∆H < 0 and with probability exp{−β∆H/λ 2 } otherwise. This is the regular Metropolis step.
The process stops when the percentage of red particles (say p 0 ) is lower than the stop parameter defined above. We point out here an additional potentially important role of the parameter λ that defines the size of the cells (block spins, cf. the terminology from [25]) in the lattice, namely λ enters as tuning parameter in the probability exp{−β∆H/λ 2 } computed in terms of β and ∆H. Essentially, we use the mesoscopic size of the cells to normalize the difference of energy with respect to the number of particles in a single cell, that is λ 2 .
Note that a Markov chain is said to be ergodic if it is possible to reach in a finite number of steps every configuration of the system starting off from an arbitrary initial configuration. Having this definition of ergodicity in mind, our system is not ergodic due to the the way we implement the evaporation part of the system. For our purposes, we assume that we can reach the configuration σ T f in a finite number of steps, hence in a finite time interval.
The mesoscopic model proposed here is very much inspired by the one developed in [15] and is an extension of the lattice model proposed by Cirillo et al. in [4]. In our earlier work, we referred to it as the λ-model. In the next sections, we explore by means of numerical simulations the capacity of our model to produce morphologies when both short distance and large distance interactions interplay. The hope is to spot genuine mesoscopic effects by comparing our simulation output with previous work done in [4,15,31] with microscopic versions of this model (i.e., when λ = 1). It is worth mentioning already at this stage that our mesocopic model can be linked with what is observed based on purely microscopic descriptions, but is unable to reach mean-field information, i.e. when λ = min{L 1 , L 2 }. One way to facilitate the latter connection, is to rescale suitably in terms of λ the structure of the Hamiltonian, including changing the factor C into C(λ) as well as considering the simulation box as Λ(λ). Inspiration from the Lebowitz-Penrose scaling indicated in [25] can be useful in our context. To understand this connection, a suitable asymptotic analysis needs to be performed. This is out of the scope of this paper.

Basic simulation results
In this section, we investigate the effects produced by the different parameters on the shape and size of the formed morphologies. The parameters involved in the simulations reported in this section are the ones independent of the interaction length scales in the system. Particularly, we analyze: • whether the phase separation is occurring or not, depending on the temperature of the system, by changing β; • how the shape of formed morphologies and the evaporation time are altered by the volatility parameter φ of the solvent; • effects of the interaction parameter between the non evaporating species, i.e. J +1,−1 , on both the phase separation and the shape of morphology formations; • how the initial ratio of species, given by the initial probabilities p −1 , p 0 and p +1 , affects the shape of morphology formations.
To perform the simulations, we build a code in Python [14,16]. We mainly use the module NumPy 4 to build and modify the lattice and the module Matplotlib 5 to print out the lattice as a figure. In order to handle a large amount of simulations that require an intensive computational effort, we run our code on the supercomputer Kebnekaise, provided by High Performance Computing Center North (HPC2N 6 ) of the Swedish National Infrastructure for Computing (SNIC 7 ).
Performing the simulations allows us to reach different morphology shapes at different iterations capturing the large time stationary behavior of our system. By displaying graphically the obtained morphologies, we want to understand the role our parameters play in the process of evaporation of the solvent through the surface of a thin film (described by our lattice), corresponding to the vertical cross section of a 3D box. We choose the reference set of parameters to correspond to what is investigated in [31] by means of a microscopic lattice-based model involving a short range interaction energy, that is λ = 1. This allows us to test the code and facilitate eventual future comparisons of results.
In this section, we consider a lattice containing 128 × 128 particles divided in yellow (with spin -1), red (with spin 0, the solvent that is evaporating), and blue (with spin +1), endowed with periodic boundary conditions. The size of the cells is fixed to be λ = 4. The initial configurations are randomly generated according to the initial probabilities p 0 = 0.4, p +1 = 0.3 and p −1 = 0.3. This choice defines an initial mixture composition of 40:30:30 that is kept in the simulations. In this case, the initial configurations are similar to the one shown in Figure 3. The tuning parameter C entering as factor in the structure of the Hamiltonian (1) is fixed to be C = 1. Moreover, the energy parameters are set to be as indicated in the following symmetric tensor: In this case, the volatility parameter is fixed as φ = 0.6 to study the effect of the temperature. In the columns we have 75%, 50%, 25% and 10% of remaining solvent from the initial amount. We notice that for high temperature (e.g., β = 0.1) we do not reach morphology formation. The required number of steps to reach the 10% of remaining solvent is 8.58 · 10 6 for β = 0.1, 8.30 · 10 6 for β = 0.2, and 8.51 · 10 6 for β = 0.6.
For the first round of simulations, we fix the volatility parameter to φ = 0.6 and look for morphology formations as β changes. The other parameters are fixed as above. In each row of Figure 4, we display snapshots of the evolution of our systems for different temperatures when the percentage of remaining solvent is the 75%, 50%, 25% and the 10% of the initial content of solvent. In this case, with a lattice containing 128 × 128 particles of which 40% is solvent, those percentages corresponds to 30%, 20%, 10% and 4% of the total amount of particles, respectively. We notice that the temperature plays an important role in the phase separation of the components: if the temperature is extremely high, i.e. in the limiting case in which β vanishes, we do not reach any morphology formation. Indeed, in this particular case the thermal energy k B T = β −1 is much larger than the spin-spin interaction term J σxσy , in the Hamiltonian (1), hence spins randomly exchange in an uncorrelated fashion.
As shown in Figure 4, we notice the first hints of phase separation if we consider a lower temperature, i.e. a higher value of β. Already for β = 0.1 we see a small decrease of the energy in the system during the evaporation process: particles are arranged in interpenetrated formations, but all the species are still mixed in those stains and sharp interfaces between the regions are not very clear/smooth. Considering β = 0.2, we notice a clearer distinction in the phase separation between blue and yellow particles. However, solvent particles are mixed inside the produced morphologies and it is quite likely that they will move away after a sufficiently long time has elapsed. Finally, if we consider a much lower temperature β = 0.6, the energy finds easier the desired minima. This results in the tendency of the solvent (red) to play the role of perimeter for the blue and yellow colored areas. Since the number of steps is fluctuating, we are unable to notice a clear trend in the number of steps needed to reach the stop parameter. In this case, we need a range from 8.58 · 10 6 (for β = 0.1) to 8.30 · 10 6 (for β = 0.8 and β = 0.2) in the number of iterations to reach the 10% of remaining solvent. The formed lanes shown in Figure 4 (last row) can be perceived as steady state, since their thickness does not change anymore. Having in mind the meaning of the volatility parameter φ, we can speed up the evaporation process by increasing the volatility of the solvent. We see in Figure 5 the following effect: as we increase φ for a fixed β, we notice a strict monotonically decreasing trend in the required number of iterations to reach the stop parameter of 10% solvent. The effect is huge already for a small φ > 0. The number of iterations goes from 827.24 · 10 6 for φ = 0, to 4.96 · 10 6 for φ = 1. We could use a polynomial of order 10 to interpolate the number of iterations but we notice that we have a better fit with the exponential function 8.17 · 10 8 · e −31.46·x + 9.70 · 10 6 (shown in Figure 5). This effect is related to the volatility step of the algorithm: if we consider a high value of volatility parameter, then the probability of the solvent to reach the top of the lattice is higher, otherwise it mostly depends on the Metropolis step, hence on the Hamiltonian (1). The role of the volatility parameter is not just to speed up (or slow down) the evaporation time. It also plays an important role for the final shape of morphology formations. In Figure  6, we study the different morphology shapes obtained for different values of φ if we fix β = 0.8. Also in this case we show (from left to right) the evolution of the system when we reach the 75%, 50%, 25% Figure 6. The temperature is fixed such that β = 0.8. In the columns we have 75%, 50%, 25%, and 10% of remaining solvent from the initial amount. We notice that for a small value of the volatility parameter (φ = 0) we get a stain-shaped morphologies, while increasing φ (already from φ = 0.1) we get vertical stripes.
and 10% of remaining solvent, respectively. In the first row we consider φ = 0. The complete absence of volatility yields the formation of a bi-continuous morphology with no preferential orientation, since the solvent is not forced to go upwards, but only moves upwards in a diffusive mode. The only way for the solvent to reach the top of the lattice depends on λ, β and the Hamiltonian (1) in the Metropolis step. This also leads to a longer evaporation time, as shown in Figure 5. If we choose φ = 0.1, we already see the effect of the volatility, resulting in almost vertical stripes for the morphologies. While for the 75% of remaining solvent we still have some stains, with the evaporation process those stains are deformed in the direction of evaporation. If we increase the volatility to φ = 0.5, we notice that the vertical stripes are thinner. A similar situation arises for the case φ = 0: in the absence of solvent, the two remaining species are not forced anymore to follow the vertical evaporation and they are reorganized following the Metropolis step. The same effect is better emphasized in the lower half of the lattice, for the case φ = 1, due to a faster evaporation of the solvent. In this limiting case, the vertical stripes are thinner, hence we obtain morphologies with predominant stripe-shaped patterns 8 . Now, we want to study which effects can be obtained when varying the value of the interaction parameter between the non evaporating species, namely when changing J +1,−1 . In Figure 7, we fix β = 0.6, φ = 0.6, the interaction tensor as and we display the evolution of the system for 75%, 50%, 25% and the 10% of remaining solvent, respectively. Similarly to the limit case when β is going to zero, we see that if we consider J +1,−1 = 0.1 J 0,+1 = J 0,−1 = 1, then we do not reach the separation of phases. In this case, the interaction between the blue and yellow particles is too weak compared to the interaction with the solvent to be able to lead to any morphology formation. Increasing this interaction to J +1,−1 = 0.9, as shown in the second row of Figure 7, we do not notice a distinct phase separation. The presence of blue particles in each yellow area and of yellow particles in each blue area is not negligible, so we do not reach a clearly shaped morphology. Looking at the case J +1,−1 = 2, we already see a strong phase separation, mainly composed of eight blue and yellow vertical stripes.
One may wonder why the morphologies are slanted in this particular case (see third row), and whether this effect is connected to some suitable combinations of model parameters. This is not the case and the observed effect is not robust with respect to changes in parameters. Such skewed orientations seem to appear due to a twofold reason: the blue and yellow particles are constrained to satisfy periodic boundary conditions and the total phase separation takes place rather fast reaching stationarity. Consequently, local agglomerations of blue/yellow particles at the top or bottom boundaries are likely to nucleate fast growing elongated morphologies.
Apparently, there is something else that plays a stronger role in this skewed formation than periodic boundary conditions: the evaporation process has a considerable dependence on initial conditions. To clarify this aspect, we show in Figure 8 the evolution for 75%, 50%, 25%, and 10% of remaining solvent for different boundary conditions on the lattice. The parameters are fixed as in the third row of Figure 7, i.e. fixing J +1,−1 = 2. In the first row we present the case with fully periodic boundary conditions, as defined in Section 2. Here, the parameters and boundary conditions are the same as the third row of Figure 7, however the slightly different initial configuration results in a completely different domain. Somehow, the main direction of the morphologies is vertical even if there is still a skewed branch connecting through the bottom-top periodicity. The last row of Figure 8 shows the evolution when bottom-top reflecting boundary conditions are used. We still consider a left-right periodicity, but we do not consider periodic boundary conditions over the top and bottom boundaries, i.e. we do not consider two sites x, y in cells X = (l 1 , x 2 ), Y = (1, x 2 ) as nearest neighbors and we do not compute the energy over the vertical periodicity if the two cells involved in the bond are close to the top or to the bottom.   close to the top boundary. As follow-up question we could explore in which way the thickness of such crust depends on key parameters like volatility or temperature.
Considering an even stronger interaction, namely J +1,−1 = 15, we get a similar effect as for high values of φ: the vertical stripes are thinner in the upper half of the lattice, while in the lower half we notice some horizontal link among the formations, since keeping more areas linked together minimizes the energy of the system. For bigger values of J +1,−1 , we see that the vertical stripes are a little bit thinner but the general structure remains the same as the case J +1,−1 = 15.
It is worthwhile to also observe the rather thick lanes of solvent particles pointed out in the first two panels of Figure 7 (first row). They appear because the amount of solvent particles in the system is sufficiently high and their interaction strength with respect to the neighboring environment (i.e. the blue and yellow particles) is high. This situation is likely not to happen in the context of organic solar cells. Interestingly though, a very similar situation like the one mentioned here occurs when large crowds of charged colloids [7] or of self-driven particles like human crowds [3] or large fish communities [9] anticipate the position of the wanted exit and collectively adapt their dynamics to reach it. The top interface, where the solvent evaporates, plays the role of the exit in such context and the red particles would be the active agents. This type of results indicate that our model is likely to find applications in the context of socioand econophysics. We will exploit alike connections with population dynamics elsewhere.
As a last simulation in this section, we study how the geometry of the morphologies is affected by choosing different initial ratios of the species involved in the mixture. We fix the temperature and volatility parameters to be β = 0.6, φ = 0.6, while the interaction tensor is as stated in (2). In Figure 9, we show the initial configuration and the corresponding evolution with 75%, 50%, 25%, and 10% of remaining solvent. We recall that the initial ratio is defined by the percentages p 0 : p +1 : p −1 , the considered cases are 20:40:40, 40:30:30 and 80:10:10. We notice that the phase separation is not affected by the initial ratio of species, while the shape of morphologies is altered. However, the evaporation time is higher if we have more solvent in the initial lattice. For a high initial percentage of solvent, we also observe that we have less links among the vertical morphologies. Starting from the case 20:40:40, we see that the formed morphologies do not follow just the vertical direction of evaporation: in this case the percentage of solvent is too low to force the other two species to align upwards, resulting in a couple of horizontal morphologies. Here, the morphology obtained for the 10% of remaining solvent is quite similar to the one with the 75%, due to a fast evaporation process. For this ratio, the 75% of remaining solvent already corresponds to the 15% of total particles. Already from the case 40:30:30, the stripes are thicker and less linked among them. In this case, the scenario at the 75% of remaining solvent is different than the end, when we have the 10%, since the evaporation is much longer than the previous case. Indeed, as specified beforehand, for the case with ratio 20:40:40 the percentage of solvent with respect to the total particles is not enough to drive the formed domain to a considerable different shape. When we consider 80:10:10 as initial ratio, we observe just eight blue and yellow main stripes in the morphology. This effect is a consequence of the high initial percentage of solvent: even if the evaporation process is longer, the upwards movement of this amount of solvent drives the morphologies to be shaped in a few vertical straight stripes. For this initial ratio, the 75% of remaining solvent is still too much to appreciate a morphology formation, since that corresponds . Temperature and volatility parameters are fixed such that β = 0.6 and φ = 0.6. The interaction tensor is the one defined in (2). In the columns, we plot the initial configuration (100% of the initial amount of solvent) and the corresponding evolution with 75%, 50%, 25%, and 10% of remaining solvent from the initial amount. In the first row, we have 20:40:40 as initial ratio of species, for the second 40:30:30, while for the third we take 80:10:10. We notice that for the third case, the vertical stripes are wider due to the considerable migration of solvent.
to the 60% of total particles (with respect to the 30% of total particles for the case 40: 30:30).
Unlike to what was seen in [4], the formed domains look rather similar at top and bottom interfaces due to the boundary conditions. Moreover, as we notice later on in Section 4, the width of the formed morphologies is clearly influenced by the size of the chosen mesoscale λ.

Multiscale effects on morphology formation
In this section, we explore the effects of varying the range of long distance interactions which can be seen at the scale of morphology formation. The parameters that we analyze are the ones involving the size of the lattice (simulation box) and of the cells (spin blocks). We probe different interaction scales for the same choice of reference parameters. In particular, we examine: • repercussions on the scale of the morphologies of the lattice size, by changing L 1 and L 2 ; • how the mesoscale length λ, defining the size of the interaction cells, affects the formation of morphology for a fixed size of the lattice; • interplay between the parameter λ and the lattice size; • effects of the tuning parameter C, arising in the Hamiltonian (1), on both the phase separation and shape of the morphologies.
These numerical experiments are meant to help us understand how to rescale the system so that increasingly larger values of λ can be taken without increasing too much the size of the simulation box. The simulations will indicate that there is a clear connection between these two parameters. Similarly as in Section 3, the results are explained by displaying the graphical output of the listed simulations and by noticing the corresponding number of iterations.
To run the simulations proposed in this section, we fix a target morphology formation, hence we make a selection of the key parameters pointed out in Section 3. Particularly, we choose the volatility φ = 0.6, the temperature β = 0.6, and the interaction tensor as defined in (2). For the initial ratio of the mixture components, we set 40:30:30. Moreover, we take as simulation box size 256 × 256 and λ = 4 as cell size. The tuning parameter entering the Hamiltonian (1) is C = 1, unless otherwise specified. In this section, we show the initial configuration in every presented figure.
Firstly, we consider the effects of different box sizes. The parameters are fixed as above, except for L 1 and L 2 . In Figure 10, we display the initial configuration (with 100% of solvent) and the evolution with 50% and 25% of remaining solvent. In the first row we consider the box size to be 128 × 128 (with 6519 red particles), in the second 256 × 256 (with 26481 red particles), while 512 × 512 (with 104881 red particles) in the third. Clearly, if we consider a lattice with more solvent, the evaporation process will be longer and will exceed a reasonable simulation length; see Section 3. Considering that in the case 512 × 512 the lattice contains an excessive amount of solvent and it took almost six days on the supercomputer to reach the 25% of remaining solvent, we avoid for this round of simulations to display the 10% column.
It may seem that if we consider a bigger box size we are just doing a "zoom out" of our system, but we also notice that the structure of the formed morphologies is slightly different. We start from the first row of Figure 10, that displays the usual morphology visible in the last row of Figure 4 and in the second row of Figure 9. Moving the attention to the second row, i.e. the case 256 × 256, we notice that the average width (in term of lattice sites) of the morphology is the same. However, in the lower half of the lattice we observe more horizontal links among the morphologies than in the upper half. As in the case shown in   the last row of Figure 9, the stripes seem to be steeper in the top half of the lattice due to the upwards movement of the great amount of solvent. Lastly, we see the same behavior if we move from the second to the last row, i.e. the 512 × 512 case. Here, the average width of the formed morphologies is still the same as in the first row. Most of the horizontal links among the vertical stripes are concentrated in the lower quarter of the lattice, while in the remaining part of the lattice the shape of morphologies is similar to the upper half of the case 256 × 256. In Figure 11 we display this sort of self-similarity type behavior: if we display the 512 × 512 lattice with the 25% of remaining solvent from the last row of Figure 10, one of its lower quadrants is a 256 × 256 lattice with a similar morphology as the last evolution of the 256 × 256 configuration in the second row of Figure 10. We can follow the same reasoning as we shift from one of the lower quadrant of this 256 × 256 lattice, that is a 128 × 128 lattice, to the proper evolution of the lattice in the first row of Figure 10. It is worthwhile to observe that this self-similarity effect holds only if the percentage of remaining solvent is small enough. A similar feature has been pointed out in Figure 10.
If we consider the evolution at the 50% of the studied cases, most of the solvent is still in the central belt of each of the presented lattices. This feature destroys similarity.
We can now switch to considering the effects of different spin block sizes, namely when altering the mesoscopic parameter λ. We chose the parameters as specified above in the section. Moreover, from the latter paragraph, we know that we can consider a bigger lattice, i.e., with box size 256 × 256, with minimal loss of generality. The size choice for this simulation is forced by the values of λ that we want to analyze. In Figure 12, we show the case λ = 2 i , i = 0, 1, 2, 3, while in the column we display the initial configuration and the corresponding evolution with 75%, 50%, 25%, and 10% of remaining solvent. Since we need λ min{L 1 , L 2 }, we prefer to use a bigger box size when it comes to the case λ = 8. We observe that the size λ can also speed up (or slow down) the process: we start with 2.79 · 10 8 iterations for λ = 1, then 2.22 · 10 8 for λ = 2, subsequently 9.73 · 10 7 for λ = 4, and lastly 3.44 · 10 7 iterations for λ = 8. For a bigger value of λ, the evaporation is faster because the particles in the lower area of the lattice can freely move from a box to another, whilst if we consider the case λ = 1, a particle in the bottom of the lattice has to "travel" across every single site before evaporating. The case λ = 1, shown in the first row of Figure 12, is widely studied with a microscopic lattice-based model and a short interaction Hamiltonian in [31]. Moving to the second row, we have the case λ = 2 that presents a thicker morphology than the previous case. We also notice less independent vertical stripes. In the third row, we can find the case λ = 4. For this spin block size, the morphologies are even wider and with more horizontal agitation. Lastly, we have λ = 8 in the last row. This case seems similar to the cases studied in Section 3, as well as the first row of Figure 10. These similarities come from a wide morphology with similar horizontal links. As we can understand from the foregoing comments, the spin block size plays an important role in the rescaling of the system, as it defines the lenght scale of interaction in the Hamiltonian (1). Particularly, increasing this size, we are decreasing the number of stripes in the formations. In fact, we are increasing the relative thickness of the morphologies with respect to the box size. It is noteworthy that if we use the spin blocks as measurement units, then the average width of the morphologies is the same for the presented values of λ.
In view of the last two paragraphs, we can spot some interconnections between the rescaling with the box size and the one with the spin block size. If we analyze Figure 10 and Figure 12 together, we can conclude that if we consider a 256 × 256 lattice and change λ accordingly, then we are able to reach morphologies that are similar to those of different box sizes. Specifically, visible correlations can be pointed out between: • λ = 4 for L k = 128, k = 1, 2 and λ = 8 for L k = 256, k = 1, 2; • λ = 4 for L k = 512, k = 1, 2 and λ = 2 for L k = 256, k = 1, 2.
We expect that, while keeping those different proportions between the box and the spin block size, we can rescale our overall system so that the structure of the morphology shapes is preserved. This would bridge information between two distinct mesoscales.
As last simulation, we examine the effect of the tuning parameter C, stemming from the Hamiltonian (1). The purpose of this dimensionless parameter is to weight the intra-cell energy with respect to the total energy. We fix the set of reference parameters as stated in the beginning of this section. We consider the cases C = 0, C = 0.1, C = 1, and C = 10. In Figure 13, we illustrate the initial configuration and the corresponding evolution with 75%, 50%, 25%, and 10% of remaining solvent. In the first row of Figure  13, we look at the results for C = 0. Here the morphology formation is striped but, since the intra-cell energy is completely neglected, we can find more solvent than usual in the morphologies, also with some impurities, i.e. blue particles in the yellow zones and vice versa. In the second row we consider C = 0.1.
This case leads to sharper formations. We also see that now the vertical stripes are really straight. In the third row, we display the case used for the other simulations when C = 1. Using this value for the tuning parameter, morphologies become visible already from the 75% of remaining solvent. Finally, the effect of C = 10, shown in the last row, is remarkable. Although the phase separation takes place, morphologies seem to meet difficulties to form coherent structures. Most of the spin blocks are filled by particles of the same species, while the location of spin blocks containing the same species is not regular enough to define a morphology. Nonetheless, we did expect this effect to happen as setting C = 10 makes the intra-cell energy disproportionate with respect to the inter-cell interaction. If we consider this tuning value in the Metropolis step, hence in the Hamiltonian (1), the inner energy of every spin block is playing a pivotal role, while the interfacial energy between two cells is almost neglected. We still observe a few vertical links among the spin blocks, but it is mainly because of the upwards movement of the solvent. Even if this situation was foreseeable, it is useful to see the impact on the overall morphology. It would be interesting to unveil a physical interpretation for C, of particular importance would be to which extent does C hold information on an eventual λ-dependence.

Discussion and outlook
The simulation tests and the corresponding discussions of the observed effects reported in Section 3 and Section 4 do not cover all possible scenarios. An exhaustive discussion of the case λ = 1 is done in [4]. Within the framework of this paper, we aimed: 1. to show the capability of the model to produce coherent morphologies at any characteristic mesoscopic length λ sufficiently smaller than the simulation box size, and 2. to explore eventual connections between simulations and morphologies obtained using different values of λ and eventually also different volumes of simulation boxes.
Our study opens a number of paths for possible further research exploiting further this multiscale model. We mention here only a few ideas which we deem as being more prominent: If one has in mind the applicability of this model to practical questions concerning organic solar cells, then one major difficulty is to set up a computable observable that can be investigated de facto in an experiment producing morphologies. Usually, as pointed out for instance in [22], in the laboratory one has access to a top view height scan of the morphologies as measured by Atomic Force Microscopy (AFM), while our simulations deliver a transversal view on these internal structures. We believe that using a kinetic Monte Carlo approach would allow a better understanding of the physical clock of the morphology formation. Moreover, extending our implementation of our mesoscopic model to a 3D version would allow a better insight as well as a good match to the AFM images, in which variations in the height of the top surface of the film are shown. For this to happen, besides extending the model, a massive computational effort is needed. Should this be successful, then the triad of methodologies (theoretical, experimental and computational) would then be integrated to make possible an adequate attack of the central two-fold question: Which morphology is best suited for organic solar cells and how to control it?
One of the main aspects that we still wish to investigate further in this context is what type of continuum models are corresponding to the mesoscopic lattice model formulated here. Particularly, we would like to explore under which conditions we can bridge suitable averages of our simulation output to what one would obtain with approximating numerically at the continuum macroscopic scale a Cahn-Hilliard-Cook-type model for a ternary mixture with the evaporation of one component. As first step, we will be investigating which of the geometric structures of the morphologies formed via our mesoscopic simulations can be obtained via changing parameters in the Cahn-Hilliard-Cook model from [11] (or other variants as reported e.g. in [28,29,32], and, more recently, in [26]), and which are not obtainable via such macroscopic-level simulations.
Once morphologies obtainable via both the λ-model as well as by the Cahn-Hilliard-Cook-type system for a ternary mixture with evaporation are classified, then one can think of studying the effect of the obtained morphologies shapes on the efficiency of the macroscopic flux responsible for charge transport. This upscaled information can be reached via averaging the transport of charges over an array of microstructures (REV) which all have as inclusions the obtained morphologies. This discussion can potentially be done at the level of the Nernst-Planck-Poisson system as in [27] and eventually it can be combined in a shape optimization framework. A similar work program has been proposed in [8], but they did not consider physically realistic selections of morphologies leaving thus place for a number of improvements.