Quantitative analysis of phase formation and growth in ternary mixtures upon evaporation of one component

We perform a quantitative analysis of Monte Carlo simulations results of phase separation in ternary blends upon evaporation of one component. Specifically, we calculate the average domain size and plot it as a function of simulation time to compute the exponent of the obtained power-law. We compare and discuss results obtained by two different methods, for three different models: 2D binary-state model (Ising model), 2D ternary-state model with and without evaporation. For the ternary-state models, we study additionally the dependence of the domain growth on concentration, temperature and initial composition. We reproduce the expected 1/3 exponent for the Ising model, while for the ternary-state model without evaporation and for the one with evaporation we obtain lower values of the exponent. It turns out that phase separation patterns that can form in this type of systems are complex. The obtained quantitative results give valuable insights towards devising computable theoretical estimations of size effects on morphologies as they occur in the context of organic solar cells.


I. INTRODUCTION
Morphology formation is one of the key factors in the processing of multi-component thin films from solution. In applications such as organic solar cells the photoactive layer is a thin film of electron donor and electron acceptor molecules. The internal structure of the morphology of this layer plays a very important role for the electric charge generation and collection; cf. e.g. [41]. The photoactive layer is produced within a solution with organic solvent -a distinguishing feature of the organic solar cells compared to the more conventional, siliconbased, photovoltaic systems. Morphologies are formed by phase separation of the electron donor and electron acceptor molecules during the evaporation of the solvent.
In the framework of this paper, lattice spin systems are used to understand size effects on different morphologies as well as characteristic time scales observed in domain growth phenomena. A paradigmatic model is the widely studied 2D Ising model [30], that served as a model system for the investigation of phase separation in magnetic materials and of novel transport mechanisms under nonequilibrium conditions [13]. For such a model, it has been well established that in the spinodal decomposition regime, namely for zero external field and subcritical temperature, the average diameter of the formed domains scales with time as a power law, where the exponents depend on the details of the dynamics. For nonconservative and conservative dynamics the exponents are 1/2, and respectively, 1/3; see e.g. [5]. Interestingly, these results were shown to be robust with respect to slight modifications of the Hamiltonian describing the equilibrium properties of the system; we refer the reader for instance to [7]. On the other hand, different interesting phenomena are observed when the dynamics are modified more drastically by introducing, for example, shear effects as in [8].
In order to describe the formation of internal structures in the presence of an evaporating solvent as typical to applications in the context of organic solar cells, a three-state model is needed. Two straightforward generalisations of the Ising model with nearest neighbour interaction between three-state spin variables are the Blume-Capel [3,4,6] and the Potts [31] models. Both these models have received much attention for their ability to model different physical situations both at equilibrium [40] and out of equilibrium [9][10][11]20]. The distinguishing feature between these two models is that in the Blume-Capel model, interfaces between different spins have different costs. The dynamics of phase separation for the Potts model seen in the spinodal decomposition regime are not completely understood: a rather clear scenario is described for nonconservative dynamics for three and four state spins [37], whereas the understanding is only partial when spin variables with larger cardinality are considered [23]. To the authors's knowledge, no references are available for the Potts nor the Blume-Capel model with conservative dynamics.
Three-state lattice models were used in [12] and [34] as an efficient tool to study phase separation patterns in ternary mixtures that allow the evaporation of one of the components as an alternative to phase-field models used in the same context [27,32,33,35]. In our previous work [12], we focused our attention on exploring the relevant parameters that influence the morphology formation. This is a subject of large interest in the community and was discussed in many experimental and computational studies [1,2,15,18,21,22,25,27,28,32,35,36]; see also [17] for related work done for stochastic models for competitive growth of phases. In the present work, the emphasis falls onto the quantitative analysis of the results for relevant choices of parameters. Quantitative studies are important particularly when one wants to compare morphologies obtained with different computational models or to compare computational and experimental results. In general, different models and computational methods give results on different length scales. In some approaches the characteristic length scales are incorporated into the model parameters, such as the parameters arising in the structure of the interaction potential for atomistic molecular simulations. Even in apparently scale-free models, the nature of the modeling assumptions induces a length-scale range for the results. Having in view potential applications of our three-state lattice models, it is of a primary concern to identify relevant quantities that capture the essence of domain size evolution across the scales. The quantitative analysis presented here is an important step towards a more general quantification and eventual classification of morphology pictures obtained by different experimental and computational methods.
In this paper, we follow up the ideas developed in our previous work [12] and give a quantitative study of the phase separated domains using a slightly simplified version of our original model generalising the Blume-Capel and Potts models. The three possible states of the spin variable are denoted here by −1, 0, and +1: "0" is interpreted as an interaction site of solvent molecules, whereas "±1" represent interaction sites of the other two components. We prefer this more general formulation, since we do not take in account the different molecular weights or volumes of the three components and hence each interaction site occupies the same volume in the lattice. In case the molecular weights of the three components are comparable, the spin variables could be interpreted as different molecules, but in our case of interest, where the non-evaporating molecules, e.g. polymers, have different sizes and are much larger than the solvent molecules, this assumption does not hold.
To study the phase separation in the ternary mixture upon evaporation of one of the components, we consider a spin model with a Kawasaki-like dynamics [24] governed by the Metropolis algorithm [26] to account for energy differences associated to possible spin exchanges and computed using appropriate boundary conditions. The Kawasaki dynamics is modified here such that evaporation of the 0 component is allowed. Keeping track of the solvent evaporation is crucial for the study of the morphology formation in solution-borne thin films, used e.g. in the preparation of the active layer in organic photovoltaics. To implement such a mechanism, the 0s in the first row of the lattice are removed (evaporation) and replaced by a +1 or a −1 with probability chosen proportionally to the initial fractions. The dynamics start with a randomly chosen configuration consisting of a fixed fraction of the three different spin species and evolve until an a priori small concentration of solvent is reached. The problem we study is related to the spinodal decomposition with the new ingredient of the evaporation of the zero component. We have to keep in mind that our model relies on assumptions and will not correspond to all aspects of the physical reality and that further improvements are possible to make the model more realistic. This falls outside the scope of the present work.
The paper is organised as follows: we first define the model and describe the tools that we rely on to investigate the domain formation under evaporation. Then, we discuss our numerical results firstly in the absence of the evaporation, and then, when the evaporation of the solvent, i.e. of the 0 component, is involved. After discussing the effect of the temperature parameter on the overall system dynamics, a short summary of our findings concludes the paper.

II. MODEL AND METHODS
In this section, we firstly define the model and then illustrate the tools that we use to measure the domain growth. In the last part of this section, we briefly illustrate these methods by discussing the standard Ising case. For more details on the Monte Carlo method and its variations, we refer the reader e.g. to the monographs [14,29].

A. Model
Let Λ be the square {1, . . . , L} 2 endowed with periodic boundary conditions. An element of Λ is called site and two sites are said to be nearest neighbours if their Euclidean distance is one. A pair of nearest neighbouring sites is called a bond. We associate the spin variable σ(x, y) ∈ {−1, 0, +1} with each site (x, y) ∈ Λ and define the total energy using the Hamiltonian for any configuration σ ∈ {−1, 0, +1} Λ , where the rows and columns of the interaction matrix J refer respectively to −1, 0, and +1. In this case, there is no cost associated with self-interaction (the main diagonal of J), a relatively small cost between ±1 and 0 sites, and a relatively large cost of an interface between −1 and +1. Such a choice of interaction matrix corresponds to a wellstudied parametrisation of the Blume-Capel model with zero magnetic field and zero chemical potential; see e.g. [9,11,20]. Note that this specific structure of the interaction matrix J promotes the phase separation of the components, irrespective of the presence or absence of the evaporation. Consider the integer time variable t ≥ 0. Fix the parameter β > 0 and refer to 1/β as the temperature. Fix c −1 , c 0 , c +1 ∈ [0, 1) with the constraint c −1 +c 0 +c +1 = 1. They are the corresponding fractions of −1, 0, and +1 spins in the initial configuration, that is at time t = 0.
The stochastic evolution σ t is constructed by repeating at each time t > 0 the following steps 2L 2 times: i) choose a bond at random with uniform probability; ii) if the bond is of the type ((x, L), (x, 1)) and σ(x, L) = 0, then replace the spin zero at the site (x, L) by +1 with probability c +1 /(1 − c 0 ) and by iii) otherwise, let ∆ be the difference of energy between the configuration obtained by exchanging the spins at the two sites of the bond and the actual configuration; exchange the two spins at the sites of the bond with probability 1 if ∆ < 0 and with probability exp{−β∆} if ∆ ≥ 0.
We stop the dynamics when the total number of zeroes in the system becomes smaller than L 2 /10. Based on this description, the dynamics are of Kawasaki type complemented with a Metropolis updating rule with the addition of the evaporation rule (i.e. step ii) of the algorithm).

B. Methods
A common measure of the average domain size is obtained by fixing a cutoff for the two-point correlation function. More precisely, for any s ∈ {s x , s y } let be the two-point correlation function. Moreover, we also consider the horizontal and vertical two-point cor- , 0), t) and G y (r, t) = G((0, r), t), where r is an integer number. The directiondependent two-point correlation functions typically decrease from their maximum value at s = 0 in an oscillatory fashion, such that it is possible to estimate the size of the domains by fixing a cut-off and finding the value at which the two correlation function intersect such a cut-off. In this way, we shall find an estimate for the horizontal and vertical diameters of the domains, R x and R y , respectively. Another well-established domain-size measurement is based on the first momenta of the structure factor. More precisely, for any (k x , k y ) in the first Brillouin zone {−π, −π + 2π/L, . . . , π − 2π/L, π} 2 , let be the structure factor. Note that the quantity in the absolute value above is simply the Fourier transform of the configuration at Monte Carlo time t, and can thus be evaluated using any form of fast-Fourier-transform technique to speed up execution. Hence, another way of estimating the horizontal and the vertical diameters of the domains is by redefining R x and R y as follows: where α ∈ {x, y} and each summation is carried out over the first Brillouin zone.

C. Domain growth in the 2D Ising model
As a test of this methodology, we consider the classical two-dimensional Ising model under Kawasaki dynamics, for which the growth exponent of 1/3 has been extensively verified; we refer the reader for instance to the works [5,7]. The results for our case are shown in the left and right panel of Fig. 1. Here, we report the resulting domain growth of the same run as evaluated using both the correlation function and the structure factor. The previously established 1/3 exponent is recovered, lending credit to the chosen methodology.
Based on Fig. 1, we note that the two methods yield quite similar power law structures, i.e. , the same exponent with different constant prefactors. This yields a discrepancy in the absolute values of the domain sizes calculated based on these two methods. Since we do not have a physically imposed lengthscale (i.e., no physical meaning of length is included in our lattice models), this is not a critical distinction. It is more important to observe how the domains grow, hence we wish to compare the exponents in the power laws mainly at long times. Sometimes, also short timescales could be of interest if power laws of suitable observables are detected. Both the correlation function and structure factor measure show that two different regimes can be distinguished: the initial one in which the domains start to be formed by coalescence of equal spins and the second one, characterised by the power law scaling, in which the already formed domains grow in time. This last regime will be addressed as the growing or scale invariant regime.

III. NUMERICAL RESULTS
After having checked the validity of the different methodologies on the Kawasaki Ising dynamics, we now study the three-state model introduced above. We explore the model for different choices of the initial fraction of zeroes c 0 . For all the simulations, we set c −1 = c +1 , so that the initial number of minuses and pluses will be equal. Due to the evaporation mechanism, the ratio between minuses and pluses will oscillate slightly during the overall evolution, while the fraction of zeroes will progressively decrease.

A. Domain growth in the three state model without evaporation
The results for the growth of domains obtained with the two different methods for calculating R x from the same Monte Carlo simulation of the three state model are shown in Fig. 2. We omit to show the results for R y as in the absence of evaporation domains are isotropic, as in the Ising case (see Fig. 1). Data for different values of the initial zeroes concentration c 0 are reported here. The connecting lines between the data points do not imply piecewise linear regression. The x 1/3 line is shown as reference, since a 1/3 exponent of the power law for the domain size as a function of time is the observed behaviour in the scale invariant regime via the correlation function (left panel). This is consistent with the findings for a conserved dynamics as reported in [5]. One single system history for each solvent concentration c0 is considered. The two panels show that in the two state case, namely, for c 0 = 0, growth is very slow and in the time interval considered in the simulation, domains have just started to form. This simply means that the scaling regime has not yet been reached, and that the overall process is still in its incipient phase. This line of reasoning is supported by comparing the configurations shown in the upper two rows of Fig. 3. Indeed, looking at the first row, configurations referring to the case c 0 = 0.0 are shown. They point out clearly that the growth regime has not yet started, compared to the second row with c 0 = 0.1, where the domains have grown substantially.
Both analysis techniques agree when predicting that the domain growth is much faster when a moderate amount of zeroes is present in the system. This situation is not unprecedented in literature, see, e.g., [19] where the authors demonstrate a process of greatly speeding up the Ising Kawasaki dynamics by introducing one or several vacancies in the lattice. We argue that the zero component of the three state model acts as a form of vacancy, since its interface cost is smaller. Hence, replacing it by a minus or by a plus spin is cheaper. Thus, zero sites are more mobile and greatly speed up the dynamics. Indeed, for c 0 = 0.1, which is reported in the second row of the same figure, the domain formation started at about 10 3 MCS. This is also very well confirmed by the data in the right panel of On the other hand, the structure factor-based measure suggests that growth is slower when the zeroes concentration is increased. The configurations plotted in the lower three rows of Fig. 3 indicate that the growth mechanisms change when more zeroes are present. While in the second row of Fig. 3 zeroes form a thin film around plus and minus domains, so that growth happens essentially as in a two-state lattice system, in the plots provided in the lower three rows of Fig. 3 the situation is rather different: domains of zeroes have sizes comparable with pluses and minus ones. Furthermore, the three species seem to compete during growth. In the last row, when the zero concentration is very high, the process seems to become more peculiar, in the sense that plus and minus domains grow inside a connected background of zeroes. It is worth mentioning that similar behaviors have been observed in different regimes, specifically in the study of the metastability occurring in the framework of the Blume-Capel model [9][10][11]20], when growth does not happen via the coalescence of small droplets but rather via a sudden nucleation of a large droplet.
The fact that different values of the zeroes concentration give rise to different growing mechanisms suggests that, for this three state system, the way in which the structure factor measures the size of the domains is more reliable than that provided by the correlation function technique. Even if, from the experimental point of view, this case might seem less interesting, our results can be related to situations where the concentration of the active components in the solvent does not vary much, as it would be for instance the case for a certain short amount of time at the bottom of the film or for a very slow evaporation.

B. Domain growth in the three state model with evaporation
Now, we discuss our results for the domain growth in the presence of the evaporation, which is the most interesting case in terms of applications to ternary mixtures in the context of organic solar cells referred to in the introduction. The horizontal and vertical domain sizes measured with the correlation function and with the structure factor are reported in Fig. 4 and Fig. 5, respectively. Just like in the case discussed in the previous section, the two-point correlation function shows a domain growth exponent of 1/3 in the appropriate regime, while the structure factor method gives a more complex behaviour, showing growth exponents between 1/5 and 1/3 for R x . What concerns R y , we see briefly a growth exponent of 1/3, which then speeds up further near the end of the evaporation. This asymmetry is not necessarily surprising as this problem is anisotropic due to the evaporation of the zeroes at the top-row of the lattice.
Even in the presence of evaporation, the growth mechanism seems to depend on the solvent concentration, see Fig. 6. Note that the last frame (for 729415 MCS) is missing for the first two initial concentrations of zeroes, as the length of the simulation is decided by the final concentration of zeroes (recall the stopping condition for the dynamics is when c 0 = 0.1) and, in these cases, that value is reached for shorter simulation lengths. By similar arguments as before, it is believed that the structure factor yields the most meaningful domain size calculation, which is further supported by the observation that R y > R x near the end of the evaporative process (compare Fig. 4 and Fig. 5). This also appears to be the case when looking at the final configurations in Fig. 6, and it is thus believed to be a good indicator of the validity of the domain size calculations.

C. Effect of temperature
Temperature is an important factor in the context of our lattice models as well as for the actual experimental processing of the thin film, both in the initial phase of the solvent evaporation, but also in the late stages and even for after processing, via thermal annealing. Direct comparisons of the simulation results with experiment are not feasible at this stage due to lack of temperature controlled in-situ experimental data especially at early stages but also due to the qualitative character of the temperature in our model, captured only in terms of β −1 . Nevertheless, we can see clear tendencies on the domains growth and this is a good starting point for further investigations. In this section, we investigate the effect the temperature β −1 has on the domains growth as they are formed in the context of the three-state lattice model without and with evaporation. The estimates of the domains size are done via the structure factor method.
In Fig. 7, we plot our result in absence of evaporation. Several noteworthy aspects appear. With a binary ±1 mixture, namely, c 0 = 0.0 (see the left top panel in the figure), it is clear that increasing temperature is associated with accelerating dynamics. As the concentration of zeroes is increased, this behaviour becomes more complex. More specifically, the initial trend is still compara-ble, i.e. increasing temperature can be associated with larger domains after the same time. At larger times, this is no longer the case, see the time slice ∼ 10 5 MCS in the top right panel of Fig. 7, where the β = 1.2 data crosses the β = 0.8 one. As the solvent concentration continues to increase, the effect of the temperature diminishes before the inflection point (which also occurs at earlier times, bringing evidence to the idea that the third species speeds up the dynamics). Finally, for c 0 = 0.6, the inflection point is no longer visible in the current time domain; see the right bottom panel in Fig. 7. It is instructive to look for a moment at the morphology formation exhibited in Fig. 8 for the case c 0 = 0.4. Note that for a high temperature (β = 0.8), the zero sites (depicted in red) readily penetrate the minus-and plus-filled domains (the yellow and blue regions), to the point where the interface boundaries are rather difficult to identify. When the temperature is low (e.g. for β = 1.2), the situation is significantly different, and the phases are clearly defined, with minimal interpenetration. It is difficult to say by inspection that the domains are larger in terms of area in the second case (i.e. lower row in Fig. 8), but it is clear that they are more well defined. Since the Fourier transform is quite sensitive to the sharpness of the boundaries of the domains, this will affect the results in Fig. 7 in the sense that results for high β are more accurate compared to those for low β.
We further consider the effect of temperature when the zeroes evaporate through the top-row of the lattice, see Fig. 9. These data appear consistent with the case without evaporation at least from the point of view that the initial behaviour is similar. At longer times, the evaporation effects dominate, and the domain coarsening deviates from scenarios computed with the three-state lattice model without evaporation. Here again, we observe the inflection point, but only in the first row of Fig. 9, indicating that the dynamics are further accelerated due to the evaporation.
A proposed mechanism explaining the inflection points in Fig. 7 and Fig. 9 is the following: high temperatures and/or high ratio of ±1 sites favours the growth of the domains. This seems to hold until a certain critical domain size is reached, after which both the increase in temperature, and hence, increase in fluctuations allow the zeros to penetrate the domains more easily. Consequently, this leads to an altering of the average domain size.

IV. CONCLUDING REMARKS
In this work, we proposed a quantitative analysis of phase formation and domain growth in ternary mixtures, both with conservative and non-conservative dynamics. In the latter case, one component is evaporated from the top of the lattice in a process akin to that utilised in the fabrication of solution borne thin films used in photovoltaics applications based on organic solar cells. We reproduce the 1/3 exponent for the Ising model, while for the ternary-state model without evaporation and for that with evaporation we obtain lower values of this exponent. Estimating how domains grow is a quite complex task -in our context, it is heavily influenced by the initial mixture concentration as well as by the temperature. Interestingly, note that the morphologies obtained in the present simulations can be found as well in experimentally-built thin films; see e.g. [27,38].

V. FURTHER RESEARCH
Suggestions for further research include an extension to three dimensions for improved applicability to experiments and to validate the two-dimensional results by studying a slice of the morphology. Staying in two dimensions, once could conceive further of a number of interesting avenues of research, such as the influence of changes to the interaction matrix J and extending the model to a more physically relevant one by taking into account the different molecular weights and volumes of the species. In the introduction section, we referred to the experimental situation that raised the domain size question we are addressing here. It is rather difficult to monitor the domain size growth at early stages experimentally. One attempt to capture early stages of the phase separation is to prepare the evaporating films under microgravity conditions and then simulate numerically those scenarios. Alike investigation routes will be studied elsewhere.