Role of the interplay between spinodal decomposition and crystal growth in the morphological evolution of crystalline bulk heterojunctions

The stability of organic solar cells is strongly affected by the morphology of the photoactive layers, whose separated crystalline and/or amorphous phases are kinetically quenched far from their thermodynamic equilibrium during the production process. The evolution of these structures during the lifetime of the cell remains poorly understood. In this paper, a phase-field simulation framework is proposed, handling liquid-liquid demixing and polycrystalline growth at the same time in order to investigate the evolution of crystalline immiscible binary systems. We find that initially, the nuclei trigger the spinodal decomposition, while the growing crystals quench the phase coarsening in the amorphous mixture. Conversely, the separated liquid phases guide the crystal growth along the domains of high concentration. It is also demonstrated that with a higher crystallization rate, in the final morphology, single crystals are more structured and form percolating pathways for each material with smaller lateral dimensions.


Introduction
Organic solar cells (OSC) are a promising technology in the field of photovoltaics. They can be solution-processed onto flexible substrates, opening the way to straightforward, eco-friendly and low-cost manufacturing, as well as novel applications such as flexible or semi-transparent solar modules. Even if their performance stays still far behind the performance of classical silicon solar cells, the best efficiencies of OSC on the laboratory scale has rapidly increased from 10% in 2014 [1] and 13%-14% in 2017 [2] [3] to more than 16% in 2019 [4] [5] [6] . The efficiency of large area organic solar modules also quickly increases and recently passed 12%. Nevertheless, stability is a second important challenge the organic photovoltaics (OPV) community is facing on the way to commercial products. Whereas silicon solar cells are inherently very stable and have a lifetime of 20 years and more, the performance of OSC tends to drop faster with time. The state of the art is that typical lifetimes of 2-4 years can be currently reached [7] with the best reported extrapolated lifetime being about 10 years. [8] Despite these encouraging results, the stability of OSC still needs to be improved. The instability of OSC is related to various different extrinsic and intrinsic degradation mechanisms. [[7]] [9] [10] Chemical degradation of the electrodes, interfaces or photoactive layers due to UV light or the reaction with water and oxygen are known as extrinsic degradation mechanisms and can be strongly mitigated by using appropriate but expensive encapsulation methods. [11] [12] Intrinsic degradation mechanisms are due to the temperature field, as well as visible and IR light absorption inside the solar cell and can lead to the formation of charge blocking layers at the interface or to the chemical evolution of the photoactive layer (PAL) itself. For instance, photodimerization of PCBM ( [6,6]-phenyl-C61-butyric acid methyl ester) has been shown to lead to a strong loss of performance due to the reduction of the electron mobility in polymer-fullerene solar cells. [13] Morphological degradation of the PAL might also be a significant source of intrinsic instability. In general, photoactive layers of OSC are multi-phase bulk-heterojunction structures made out of at least one donor and one acceptor material. In these structures, the excitons generated through light absorption are strongly bound and have to be separated into free charge carriers at interfaces between both materials. Since the mean free path of an exciton is typically 10nm, this length scale should also be characteristic for the size of the phases. Moreover, percolated pathways to the electrodes should be available in order to efficiently extract the charge carriers. Finally, relatively pure donor and acceptor crystal phases are desired to ensure high charge carrier mobility. As a consequence, the PAL is typically composed of crystal phases of the donor and acceptor materials and an amorphous mixed phase [14] [15] like in the well-known system made of P3HT (poly(3hexylthiophene)) and PCBM. [16] [17] [18] The importance of the bulk-heterojunction (BHJ) morphology has been identified to be crucial for the OSC performance for a long time. [19] [20] [21] [22] [23] However, these morphologies are typically partially quenched far from the thermodynamic equilibrium during the deposition process and hence not stable. Therefore, they can evolve during post-processing or operation of the solar cell. For instance, strong burn-in degradation in polymerfullerene solar cells has been attributed to demixing of donor and acceptor in the amorphous phase and aggregation of fullerenes into micron clusters [24] . Different successful strategies have been proposed to overcome this problem, such as adding a second acceptor, more compatible with the donor in order to stabilize the mixed phase [25] [26] [27] [28] , the use of non-fullerene acceptors like IDTBR (rhodanine-benzothiadiazole-coupled indacenodithiophene) [29] or to kinetically quench the mixed phase. [30] These improvements are based on remarkable efforts to unravel the mechanisms of the PAL morphology formation and stability. On the one hand, they base on thermodynamic considerations such as miscibility of solvent, donor and acceptors [[25]] [ [27]] [31] [32] and the evaluation of phase diagrams [[15]] [33] [34] [35] [36] . On the other hand, the importance of the kinetic evolution towards the thermodynamic equilibrium during the deposition process and ageing, including possible transient or stable liquid-liquid phase separation (LLPS) [37] [38] [39] [40] , has been acknowledged and recently qualitatively taken into account for stability improvement [[30]] . Nevertheless, no general coherent physical framework has been proposed to understand the BHJ morphology formation and stability taking into account at the same time thermodynamic aspects such as liquid/amorphous phase stability and crystallinity and the kinetics of the system (kinetics + thermodynamics, crystallinity + miscibility). Therefore, the understanding of stability is still very material system-dependent. Phase-field simulations can be relevant to tackle this problem: the phase-field method is a wellestablished continuum-mechanics diffuse interface simulation framework for solving the kinetic evolution of thermodynamically complex systems with many phases. [41] [42] [43] [44] On the one hand, it has been widely used to investigate the crystallization in single-material systems (metals [45] or polymers [46] ) or many-material systems (alloys [47] [48] [49] , precipitation from a solution [50] [51] ). Thereby, many crystal systems can be simulated with the multiple-field phase-field (MFPF) [[41]] [ [42]] or with the orientation-field phase field (OFPF) [ [63] . In the field of organic electronics, phasefield modelling has been used to study spinodal decomposition of the mixtures during the drying of the wet film but remains limited to a few papers. Wodo and co-workers dealt extensively with ternary systems including a polymer, a fullerene and an evaporating solvent, adding also specific interactions with the substrate [64] [65] . They gained insight into the impact of process parameters on the final amorphous two-phase structure. Michels and co-workers investigated donor/acceptor mixtures with an evaporating solvent [66] [67] [68] [69] [70] . With respect to stability, Ray, Alam and co-workers performed simulations of binary amorphous immiscible acceptor-donor systems based on the Cahn-Hilliard equation. They aimed at describing the impact of annealing on phase separation. [71] [72] [73] . Finally, both processes of LLPS and crystallization can be coupled in a very natural way within a phase-field framework. However, very few groups have been dealing with both at the same time. Zhou investigated numerically spinodal induced crystallization, showing that the phase separation could promote crystal nucleation and growth, [74] Rathi investigated the competition between crystallization and LLPS in a crystalline-amorphous binary system, [75] and Saylor and Kim simulated the evaporation of an immiscible amorphous-crystalline system in solution with application to polymer-embedded drug-release films. [76] [77] [78] In these papers however, only one of the materials is able to crystallize, and no specific interaction between crystals prevents them from merging. Analytical investigations of the interplay between crystallization and LLPS are also limited. For example, Mitra calculated the rate of nucleation for heterogeneous nucleation generated at the interfaces of a ripening phase-separated liquid mixture. [79] In this paper, we propose a new phase-field model which takes into account the miscibility of the liquid materials as well as their respective crystallization properties. As a consequence, it can handle crystallization of each material and LLPS at the same time. The impingement of single crystals is also included, so that polycrystalline structures can be investigated. Based on simulations of simple binary model systems, we illustrate how it can be used to investigate the stability of OPV photoactive layers, independent of the materials being crystalline or amorphous, miscible or immiscible in the amorphous phase. Through these examples, we outline interaction mechanisms between LLPS and crystallization, and show how the kinetic properties of the system can strongly affect the transient and final morphology of the bulk-heterojunction.

Model equations
The free energy functional A phase-field framework is used to simulate the kinetic evolution of the system towards its thermodynamic equilibrium. The system is composed of n materials, out of which ncryst can have crystal phases. We describe the system morphology with the respective volume fractions of all materials in the system, but also with an order parameter for each crystalline material, whose value varies between 0 in the amorphous phase and 1 in the crystal phase. Additionally, for each crystalline material, orientation parameter fields describe the orientation of each single crystal. For 2D simulations, one single orientation field per material is sufficient, whereby two or three angles would be necessary in 3D depending on the crystal symmetry. The 3D case is not discussed in the following, but the generalization is straightforward. Each crystal has its own orientation [π; π] and is assumed to remain constant during the simulation, while no orientation is defined in the amorphous phases. The thermodynamic properties of the system are defined with the help of the free energy functional: where V is the total volume. is the local free energy density and the non-local contribution due to the field gradients. The local part of the free energy is given by . (2) The first term on the right-hand side of the equation above represents the free energy density change upon ideal mixing, with R being the gas constant and T the temperature. Following the Flory-Huggins theory [80] , 0 is the molar volume of the lattice site considered to calculate the free energy of mixing. The molar volume of the fluid i is i = i 0 and its molar size of in terms of lattice units. The interactions between materials are represented by the second term: Matkar and Kyu proposed an extension of the Flory-Huggins theory for binary systems with crystalline materials. [81] [82] Equation (4) is simply a generalization of their theory for any number of materials. The first term is the classical Flory-Huggins interaction term where , is the interaction parameter between the amorphous phases of materials i and j. Now, the crystalline materials might have a crystal phase which interact with amorphous phases, especially in the diffuse solid-liquid interface. The second term stands for these interactions, with , representing the interaction between the amorphous phase of material j and the solid phase of material k. This term can be understood considering that can be interpreted as the proportion of material k being crystallized, so that is the quantity of solid and the amount of amorphous phase interacting with this solid. [[82]] The last term stands for the solid-solid interactions, with , being the interaction parameter between the solid phases of the materials k and j. Similar to Matkar and Kyu, we write , = √ , √ , with the coefficient c ranging from -2 for fully compatible crystals to 0 for fully incompatible crystals.
The third term on the RHS of Equation (2) stands for the free energy density of phase change, according to what is commonly used for the simulation of crystallization in metals and alloys [41] [ 42] , In the equation above, is the density of the material k and Δ , = ( , − 1) its free energy density of crystallization, whereby Lk and Tm,k are its enthalpy of melting and melting temperature, respectively. If Δ , < 0, the free energy of the crystal phase is smaller than that of the amorphous phase and the material k is prone to crystallize. There is an energy barrier in the solid-liquid phase transition when varies from 0 to 1, provided that | Note that other functional forms can be used with no considerable impact on the model behavior.
0, is the value of the order parameter for which the free energy density of phase change is minimized and can be seen as the crystallinity of the material, 0, = 1 representing a fully crystalline material. As a consequence, semi-crystalline materials can also be considered with such a model.
The fourth term on the RHS of Equation (2) is a purely numerical contribution meant to facilitate the convergence properties of the simulation: for common and physically relevant parameter sets (for instance high Nχ values for a highly immiscible amorphous binary polymer system), the expected equilibrium volume fractions in the separated phases are very close to 0 and 1, so that unrealistically small time steps have to be used for the calculation to converge. To overcome this problem, a contribution to the free energy is added if the volume fractions approach 0 and 1. We choose the following form: Hence, this numeric correction term has no impact on the properties of the system provided 1 − < < . The parameter has to be kept small in order to minimize the impact on the phase diagram of the model, especially on the volume fractions of the separated phases. Finally, the non-local part of the free energy functional describes the contributions of the concentration gradients and the solid-liquid phase change to the surface tension as: where is the surface tension parameter for the concentration gradient of material i, are the surface tension parameters for the gradient of the order parameter of material k, and 1, and 2, are the surface tension parameters for the orientation gradients of material k. The last term of Equation (8) stands for the orientation mismatch energy between different single crystals of a given material and is responsible for impingement of the crystallites provided the orientation mismatch is sufficiently large. This functional form has been proposed in previous work during the development of the OFPF model, and the reader is referred to the corresponding papers for the description of the properties of the calculated grain boundaries. [[52]] [ [53]] [ [54]] The surface tension between two amorphous phases i and j is proportional to √ + , whereas the surface tension of a solid-liquid interface also contains a contribution from the phase variation and from the orientation mismatch when two grains impinge. However, the surface tension also depends on other thermodynamic properties such as the molar volumes and interaction parameters and can be computed with standard methods described elsewhere [[42]] [ [43]] .

Kinetic equations
Since they are conserved quantities, the volume fractions obeys the celebrated Cahn-Hilliard equation, initially proposed by Cahn and Hilliard for binary mixtures [83] [84] and generalized later for multicomponent mixtures [ This is actually a set of coupled continuity equations, where the material fluxes are proportional the driving force for the system evolution, namely the gradient of the exchange chemical potential.
In Equation (9), , is the chemical potential density, defined as the functional derivative of the free energy functional: The first two terms stands for to the local chemical potential, whereas the two last contributions take into account the potential due to concentration gradients and hence to surface area variations. The symmetric Onsager mobility coefficients = have to depend not only on the diffusion coefficients but also on the local mixture composition in order to ensure the incompressibility constraint and the Gibbs-Duhem relationship. Several theories have been proposed to derive correct expressions for the flux, among which the "slow mode theory" [86] and the "fast-mode theory" [87] are the most successful. Their names come from the fact that the mutual diffusion coefficient in a binary system is controlled by the slowest component in the "slow-mode theory", while it is controlled by the fastest component in the "fast-mode theory". The controversy between both theories is not fully resolved yet. However, the fast-mode theory seems to better match experimental data and can also be derived from the general Maxwell-Stefan equations framework [88] in a consistent way. Relying on these arguments, we choose to use the fast-mode theory in this paper, which leads to the following expression of the mobility coefficients: Here, the coefficients are related to the self-diffusion coefficients , of the materials i through: 8 The self-diffusion coefficients themselves are also dependent on the mixture composition { }, but unless specified otherwise, they are kept constant in this paper for simplicity. However, diffusion processes are expected to be dramatically slower in the solid crystal phases. This is accounted for with a hyperbolic tangent based dependence of the diffusion coefficients on the order parameters: Here, , is the diffusion coefficient of the material i in the amorphous phases, is the value of the order parameter around which the mobility drop is centered, controls the intensity and the steepness of the diffusion coefficient gradient from the amorphous phase to the solid phase, and is an estimate of the overall crystallinity at a given position calculated as In practice, diffusion coefficients in the solid phases are orders of magnitudes smaller than those in the amorphous phase, so that diffusion inside the crystals is fully negligible over the whole simulation time.
The order parameters obey the classical Allen-Cahn equation, Here, Mk is the mobility coefficient for the solid-liquid interface for crystals of material k and will be called "interfacial mobility" in the following. We point out however that in the general case, the crystallization rate obtained in the simulation not only depend on Mk, but also on all the thermodynamic and diffusional properties of the system. The Cahn-Hilliard and the Allen-Cahn equations together ensure that the system progressively relaxes towards its thermodynamic equilibrium, by minimizing its free energy relative to the volume fraction and the order parameter variables.
The second term of the RHS in Equation (15) includes a contribution from the orientation mismatch (see Equation (8)), so that the evolution of the orientation fields has to be calculated as well. In classical OFPF models, the kinetics of crystal orientation is governed by additional Allen-Cahn equations applied to the orientation parameter [[52]] [ [53]] [ [54]] and both order and orientation fields are fully coupled through the last term of Equation (8). This approach has two major drawbacks: first, because of this coupling, the growth rate of an isolated crystal growing within an amorphous environment depends of its orientation, which is non-physical. Second, the interfaces for orientation parameters are much sharper than those for volume fraction and order parameter gradients, inducing a high computational cost. To overcome this drawback, we propose a much simpler and computationally more efficient heuristic procedure for the propagation of the orientation field: when a crystal is growing, the value of the order parameter in the amorphous phase around it increases; when the value of the order parameter on these surrounding nodes exceeds a given threshold value, the nodes are assumed to be crystallizing and are attributed an orientation. The orientation attributed to a given node is simply the one of the (already crystallized) neighboring node with the highest order parameter. As a consequence, the orientation of a given nucleus propagates together with its order parameter field. Within a single crystal, the orientation is thus uniform, so that | | = 0 and the orientation mismatch energy term in Equation (8) is zero. Moreover, since the orientation is undefined in the amorphous phase, the orientation gradient at a solid-liquid interface is also undefined. In order to ensure that the growth rate of a crystal in a surrounding amorphous phase is independent of its orientation, we simply set | | = 0 at the solid-liquid boundaries of the orientation fields. As a consequence, the orientation mismatch contribution is only non-zero at the boundaries between two crystals with different orientations, as desired.

General description of the simulated systems
In this paper, we investigate the time-dependent morphology of binary mixtures depending on their thermodynamic and kinetic properties. This contributes to the understanding of the stability of OSC bulk-heterojunctions: PAL dry films typically consist of a mixture of one donor and one acceptor material and have a complex topology with potentially two crystalline phases and one (mixed) or two (separated) amorphous phases, depending on the material properties. The morphology obtained at the end of the drying process, or even at the end of the whole fabrication process is in the general case out of equilibrium and evolves during the cell operation. This often leads to a noticeable loss of performance. In order to illustrate some of the various physical properties that might be encountered for different donor/acceptor combinations, we investigate three different kinds of systems: fully amorphous immiscible mixtures (as a reference system), mixtures of two crystalline materials that are miscible in the amorphous phase, and mixtures of two crystalline materials that are immiscible in the amorphous phase. We intentionally consider very simple model systems in order to highlight and focus on the physical mechanisms driving the morphology evolution and the interactions between them. Nevertheless, their material properties might be far from the ones of common experimental donor/acceptor mixtures and thus we only obtain a qualitative description of the involved physical processes. Performing quantitative simulations for OPV blends however requires careful measurement of the thermodynamic and kinetic properties of the mixture and will be the topic of future work. The parameters used in the simulations are summarized in Table 1 to Table 5 unless specified otherwise in the text. With a focus on the stability of the PAL, the starting point of the simulation should be the morphology that arises from the fabrication process. Although a global qualitative picture of this morphology has been identified over the years, the microscopic details remain in general unknown. As a consequence, the starting point for our simulation is a homogeneous 1:1 mixture (unless specified differently in the text) with an initial random perturbation of the concentration and (for crystalline systems) 20 randomly distributed nuclei with a diameter of 12nm of each material. Further additional nucleation is not taken into account. These nuclei are sufficiently large to be stable, i.e. the free energy change upon crystallization dominates the contribution of the surface tension, so that they are expected to grow. Thus, in particular for immiscible systems, the state corresponding to the PAL morphology at the end of the process is not the starting point of the simulation, but corresponds to a later moment, typically when the spinodal decomposition has already started. We assume that the PAL morphology corresponds to the structure of our simulated system shortly after the start of the simulation and that the evolution we describe mainly corresponds to the PAL evolution after the drying process.
Note that we perform simulations in the situation where Δ , < 0 and with stable, growing initial nuclei. Nevertheless, during crystal growth in a multiphase and polycrystalline mixture, interphases with high local curvatures might arise and ⁄ might be locally negative due to the surface tension term, leading to local dissolution of the crystal. However, since crystal surface growth/melting is a thermally activated process and since the free energy of the crystal phase is much smaller than the one of the amorphous phase, the rate of crystal melting is expected to be negligible as compared to the rate of crystal growth. To take this kinetic effect into account, in this work we simply set Mk=0 whenever and wherever the RHS of Equation (15) becomes negative. The equations are solved in 2 dimensions with a system size of 512x512 nodes. The mesh size has been adapted so that the thinnest encountered interfaces are at least 5 mesh points thick and fixed to 2nm unless specified differently. The equations are numerically solved using an Euler explicit finite difference scheme.

Amorphous immiscible systems
The case of LLPS phase separation in binary amorphous systems has been studied extensively theoretically [89] [73]] . In this section, we illustrate very briefly this situation in order to provide a reference situation for the case of crystalline immiscible systems described later. We present simulation results for symmetric simple small molecule model systems (see Table 1) and for strongly asymmetric polymer-small molecule systems (see Table 2). The phase diagrams of these mixtures are shown in Figure 1, whereby the simulated systems are marked with a black star. Additionally, for the polymer-small molecule system, the diffusion coefficients are assumed to be composition dependent. Several models have been proposed in the literature for the expression of these coefficients [92] . We propose here to use the equation proposed by Vignes that has both advantages of being a good first order approximation of the well-known dependence of the self-diffusion coefficients of polymers in solution, and of expressing the selfdiffusion coefficient depending on the self-diffusion coefficient at infinite dilution , →1 , which are experimentally more accessible: The diffusion coefficients at infinite dilution are set to be 10 -13 , 2 . 10 -11 , 2 . 10 -12 , 2 . 10 -10 m 2 /s for the polymer in polymer, polymer in small molecule, small molecule in polymer, small molecule in small molecule, respectively.    Table 3: equilibrium properties of the amorphous systems The thermodynamic parameters and composition of the mixtures are chosen such that they are all unstable and thus demix through spinodal decomposition (Figure 2a). The composition of the separated phases as well as their respective proportion in the demixed system can be readily obtained from the bimodal curve and the lever rule and are reported in Table 3. Remember that if the morphology presented in Figure 2 for the symmetric system has co-continuous pathways, this is not a general rule: the ability to obtain co-continuous pathways instead of isolated domains of the minority phase into the majority phase depend on the proportion of both phases in the demixed state, and hence from the blend composition, but also on the molecular size and on the interaction parameter. The generated amorphous phases then grow by Ostwald ripening (Figure 2b). At late stages, the characteristic domain size L-L0, with L0 being the size at the onset of spinodal decomposition, is known to increase in a binary system as − 0~( Λ 11 2 2 √ Δ ( − 0 )) 1 3 ⁄ , where t0 is the time at when the spinodal decomposition onsets and Δ the change in free energy between the mixed and demixed states. This is shown in Figure 3, where the characteristic length is calculated as follows: the time-dependent structure factor is calculated as the 2D-Fourier transform of the volume fraction field, and integrated over all directions to get the probability distribution ( , ) of q-vectors at each time. The characteristic length scale is the inverse of the mean value of q over this distribution, ( ) = 1 ∫ ( , ) ⁄ . Note that this scaling law also holds for the polymer-small molecule system with composition-dependent diffusion coefficients.

Crystalline miscible systems
In this section, we focus on crystalline miscible systems. Although the thermodynamic properties of the mixture are of primary importance in order to understand the behavior of real OPV systems, we highlight here the importance of purely kinetic properties and the mechanism of crystal development in the evolving structure. The objective is to analyze the interplay between the single crystals and the balance between diffusion rate and crystal growth rate. These findings are of noticeable importance for the time-evolution of classical OPV systems that are composed of three phases with a stable mixed amorphous phase, but also for the understanding of the more complex crystalline immiscible systems (see next section). The model system we study consists of two miscible materials with identical thermodynamic and kinetic properties. The model parameters are summarized in Table 4. The thermodynamic parameters are chosen so that the solid phases are strongly immiscible, and that both materials are fully crystalline and highly pure in the solid phase. However, the amorphous phase is thermodynamically stable, whatever its composition. The thermodynamic stable structure is expected to consist 100% of pure crystals from both materials. In order to analyze the influence of the competition between crystal growth and diffusion processes in the amorphous phase, the parameters are adjusted so that the crystal growth rates are sufficiently high to generate concentration gradients in the amorphous phase.   Before turning to the polycrystalline case, we investigate the growth rate of a single crystal in a mixture of both materials. In that case, the crystal growth rate is constant and it has been established for long that in a pure material, the Allen-Cahn equation results in a growth rate that is strictly proportional to the interfacial mobility. [[42]] [ [45]] Here, we calculate the crystal growth rate depending on composition for different M values, and plot the ratio of the interface velocity to M depending on the volume fraction in Figure 4. As expected, the ratio is independent of M in the pure material (φ=1). For all values of M, the crystal does not grow for volume fractions smaller than roughly 25%: below these values, the cost for the increase of the crystal surface due to surface tension is higher than the gain in volume due to the lower free energy in the solid phase, and the nucleus is not stable. Above 25%, the growth rate increases with volume fraction. There is a strong deviation from linearity for low interfacial mobility M, that is more pronounced around 50%. This is due to the comparatively stronger 1 2 1 2 12, solid-liquid interaction term for mixtures with φ close to 0.5 in Equation (4). For higher M values, the ratio of crystal growth rate to interface mobility becomes lower, except in the pure material. In this regime of high interfacial mobility relative to the diffusion rate, the crystal growth becomes diffusion-limited. This is a consequence of concentration gradients that appear around the crystal, because the diffusion process in the amorphous phase is not fast enough to compensate for the material consumption at the crystal surface. The effective concentration at the surface is therefore reduced, and so the growth rate. We now investigate the case of a polycrystalline 1:1 mixture with the parameters summarized in Table 4. The evolution of the volume fraction field for the 1 st material, superimposed with the location of the crystallites (yellow for material 1, dark blue for material 2) is shown in Figure 5. Spherical crystals develop isotropically until they impinge together (Figure 5a and b). The crystal growth is then determined by the available space between the crystals (Figure 5c to f). At some point (Figure 5d), the solid crystals have almost quenched the topology of the material 1-rich and material 2-rich zones, and the subsequent evolution of the system consists simply of the crystallization of the remaining amorphous domains. The morphology of the final structure (size, topology of the crystals) depends strongly on the location and number of the initial nuclei. The crystal growth rate at the beginning of the structure evolution, evaluated through the evolution of the mean equivalent diameter, is exactly the one expected from the single-crystal simulations discussed above (see Figure 6). However, the growth rate starts dropping as soon as the crystals impinge. A second contribution to this slow-down is that concentration gradients in amorphous areas surrounded by several crystals are higher than in the single crystal case. This example shows how not only the thermodynamic properties and the relative speed of crystal growth and diffusion processes, but also the overall composition of the mixture and the nucleus density have to be taken into account in order to predict the final film structure. Figure 6: comparison of time-dependent crystallite growth rate for a monocrystalline system and a miscible polycrystalline system (mean equivalent crystal radius), 1:1 blend, M=1e5

Interactions between crystal growth and spinodal decomposition in immiscible systems
In this section, we focus on crystalline immiscible systems, such as P3HT with PCBM, or PCE11 (PffBT4T-2OD) with PCBM, for which the amorphous mixed phase has been shown to demix generating a strong burn-in degradation of the solar cells. [[15]] [ [25]] [ [27]] The amorphous phases quickly undergo a spinodal decomposition and the separated liquid phases grow through Ostwald ripening while crystals grow at the same time. On top of the previously described processes, the LLPS and the crystal growth influence each other. The objective is to understand the interplay between them and to analyze how the crystals grow in the phase separated fluid providing strong and time-dependent concentration gradients.
Once again, the model system we study consists of two immiscible materials with identical thermodynamic and kinetic properties. The model parameters are summarized in Table 5. The thermodynamic parameters are chosen so that the liquid phase as well as the solid phases are strongly immiscible, and that both materials are fully crystalline and highly pure in the solid phase. The thermodynamically stable structure is expected to consist exclusively of pure crystals from both materials. In order to illustrate the importance of the kinetics at fixed thermodynamic properties, the crystal growth rate is varied, all other parameters being fixed. This allow us to investigate the influence of the competition between crystal growth and diffusion processes in the amorphous phase on the structure development and final film morphology. Note that to this end, the parameters are adjusted so that the crystal growth rates are relatively high, in the same order of magnitude compared to the Ostwald ripening rate.

Crystal development
The evolution of the volume fraction field for the 1 st material, superimposed with the location of the crystals (yellow for material 1, dark blue for material 2) is shown in Figure 7. The spinodal decomposition onsets at a very short time. Remarkably, the amorphous separated phases are organized in concentric zones around the crystal (Figure 7a and b), contrary to the amorphous case (compare with Figure 2). This is because the spinodal decomposition is triggered by the concentration gradients around the crystals and not by the Gaussian fluctuations in the amorphous phase. Then, since the growth rate is concentration dependent, the crystals tend to grow in the highly concentrated liquid phases, according to their geometry (Figure 7c, d and e). Thus, the crystal growth is guided by the constantly evolving topology of the LLPS. Although the interfacial growth rate is fully isotropic, crystals acquire a highly structured, possibly branched morphology depending on the pathways they find for growth. Since crystals of a given material all grow along the amorphous highly concentrated phases that usually forms percolating pathways for this 1:1 blend, the crystals belonging to a single material also tend to form percolating pathways (Figure  7e and f). However, at the same time, growing crystals constantly capture material from the amorphous phases and interfere with their ripening. They may create separations in the amorphous percolated pathways, and they finally completely quench the ripening when the crystallinity is sufficiewntly high to hinder any evolution of the amorphous phases (Figure 7e and f). The final state (Figure 7f) looks like a classical spinodal LLPS pattern, but this is a solid structure that does not evolve any further. Note again that in all the cases investigated here, the crystal growth rate is high enough so that crystals grow significantly during the Ostwald ripening of the amorphous phase. If the crystal growth rate would be very small compared to the rate of Ostwald ripening, the spinodal decomposition would start as depicted in Figure 7b and c, and then would develop very similar to the fully amorphous case (see Figure 2 and Figure 3). The time-dependent growth rate of single crystals shows also remarkable features at the beginning of the structure evolution (see Figure 8). We have already outlined that in a pure material system with a single crystal, the interface speed is constant with time and proportional to the interfacial mobility, so that the crystal radius growth linearly with time (full lines in Figure 8). In a binary polycrystalline system, three distinct growth phases can be recognized in the evolution of the mean equivalent crystal radius (dashed lines in Figure 8). At very short times, just after the onset of the spinodal decomposition, crystals absorb all the material around them so that they are quickly surrounded by a circular depleted liquid zone. The growth rate strongly drops, which is the reason for the « plateau » in the crystal radius evolution at small time. This depletion zones disappear with the ripening of the liquid phases which feed the crystal growth. This effect is less pronounced at low interfacial mobility because crystals grow too slowly to generate a depletion zone around them. In a second phase, single crystals grow isolated from each other in the highly concentrated amorphous phases. The crystal growth rate is very high and close to the growth rate in a pure monocrystalline system, because the separated phase in which the crystal grow are almost pure with the chosen thermodynamic parameters. During this phase. This holds until the first impingement between crystals. Then, the growth rate progressively decreases with increasing steric hindrance and progressive lack of remaining amorphous material to feed the crystal growth. The overall crystallization rate (defined as the derivative of the total crystalline volume with respect to time) is shown in Figure 9. Remarkably, the higher the interfacial mobility, the more irregular the overall growth rate. This is because for high interfacial mobilities, crystals are more likely to consume the material in highly concentrated amorphous phases around them, creating depletion zones and a sudden drop of their growth rate, until a new highly concentrated phase form around them and the growth can be very fast once again. At low interfacial mobility, the growth rate becomes much smoother, close to the situation of the miscible system investigated in the previous section (black curve).

Influence of the crystallization rate on the final structure
In this section, we study the impact of the interfacial mobility on the final, stable film structure. The order parameter of the final state for interfacial mobilities ranging from M=3 . 10 4 s -1 to M=10 6 s -1 are shown in Figure 10. Here, in order to visualize the order parameter of both materials on the same figure, the order parameter field for material 2 has been set to vary from 0 (amorphous) to -1 (crystalline) and both order parameter fields are added together. The higher the interfacial mobility, the more structured the final single crystals and the polycrystalline morphology: with higher interface mobility, crystals develop faster and hence following a still very finely interpenetrated amorphous phase separated mixture. In parallel, the final structure is finer because the crystal kinetically quench the system before the liquid phases have time to coarsen. This can be observed by plotting the time-dependent characteristic wavelength of the system computed from the structure factor of the volume fraction field (Figure 11). Whereas the characteristic size of the domains (L 3 -L0 3 ) 1/3 increases as t 1/3 in an amorphous system, it suddenly reaches an asymptotical value in crystalline systems, when the presence of the crystals quenches any further ripening of the amorphous phases. This quench occurs at shorter times, and therefore generates shorter characteristic sizes if the interfacial mobility is higher. The highest characteristic length is reached in the miscible system investigated in the previous section, for which the crystal growth is not influenced by the LLPS. Figure 11: Characteristic wavelength computed from the structure factor of the volume fraction field for an amorphous symmetric system, a crystalline miscible systems and crystalline immiscible systems with different interface mobility; t0 is the time for the onset of spinodal decomposition and L0 the characteristic length scale at t0 Furthermore, it can be observed that for high interfacial mobility (Figure 10a and b), some areas remain amorphous, although it would be expected from thermodynamic parameters that the whole system would crystallize. In fact, these areas can be considered as amorphous "defects" in the crystalline structure that are kinetically generated. They can form for instance if an amorphous volume happens to be fully surrounded by a crystal of the other material, thus being inaccessible for its own crystallization to take place (unless new nuclei might form, which has not been taken into account in this work). This is the case for the simulation presented in Figure 10b. They also can form even if a crystal could in principle topologically reach the remaining amorphous zone, but should go through a bottleneck (Figure 10a). This bottleneck generates a surface tension that is too high for the crystal to grow further. Even if the morphology is not in thermodynamic equilibrium, they are long-living metastable structures, the defects being kinetically quenched as long as the solid crystals can be considered as fixed. It should be emphasized that the appearance of the defects depends on the location of the initial nuclei and is not systematic. Nevertheless, the probability of existence of such defects increases with interfacial mobility. Figure 12: Crystalline volume of the 1 st material related to the whole volume for a crystalline miscible systems and crystalline immiscible systems with different interfacial mobility Figure 12 shows the crystallinity of the 1 st material (the evolution for the 2 nd material is similar), defined by the total crystalline volume of this material normalized by the total system volume. In such a 1:1 blend, the final crystallinity is expected to be a little bit less than 50%, taking into account the volume of the not fully crystalline interfaces. The crystallinity reaches its final value after the characteristic wavelength has stabilized, meaning that the phase morphology is fully quenched before the crystals are completely grown. The presence of the defects can be recognized in the lower final crystallinity values at higher interfacial mobilities. At last, Figure 13 shows the time needed to reach the final morphology, which is a good indicator of the stability performance of the mixture, depending on the interfacial mobility. It is found that it follows a power law close to M -3/4 . Further investigation and simulations on more systems with different parameters are needed to investigate the generality of this finding, and to understand how the time to equilibrium depends on the other properties of the system.

Conclusions and perspectives
The kinetic evolution of immiscible crystalline mixtures has been rarely investigated theoretically in the literature. However, it is of primary importance in order to understand the evolution of the photoactive layers of organic solar cells during post-processing and operation, so that the stability of these bulk-heterojunctions can be improved. In this paper, we presented a new phase-field simulation framework for the investigation of such systems. The free energy functional contains the description of the mixing term through a Flory-Huggins contribution, as well as the free energy change upon crystallization for each crystalline material and gradient terms generating surface tension effects. This allows to handle liquid-solid phase changes and liquid-liquid demixing at the same time. The volume fraction fields and the order parameter fields evolve towards the thermodynamic equilibrium via the Cahn-Hilliard and the Allen-Cahn equation, respectively. Additionally, the orientation of each crystal is taken into account, generating crystal impingement through a misorientation energetic contribution in the free energy. The proposed simulation code is three-dimensional, although only 2D simulations have been presented in this paper, and it can handle any number of materials. We used this model to investigate the evolution of crystalline miscible and immiscible binary systems, in order to better understand the stability of bulkheterojunction photoactive layers.
In systems where both donor and acceptor materials are crystalline but miscible in the amorphous phase, the mechanisms of polycrystalline growth have been highlighted. Compared to the growth rate of a single crystal in a pure material, the growth rate in the binary system drops with the concentration. It can also become diffusion-limited if the interfacial growth rate is high enough because depleted zones then surround the crystal, limiting the amount of material available for growth. Finally, impingement between crystals limits the growth rate and are responsible for highly anisotropic single crystals, although the growth rate is fully isotropic. In immiscible binary systems, the liquid-liquid demixing and the crystallization processes strongly interact together. Starting from a homogenous amorphous phase, the spinodal decomposition is induced by the crystal nuclei and the amorphous separated phases organize around them at short times. Material consumption through crystal growth perturbs the Ostwald ripening of the amorphous phases. With increasing crystallinity, this Ostwald ripening becomes progressively quenched by the presence of the solid crystals. Conversely, the crystal growth is induced by the spinodal decomposition: it is driven by the local concentration, so that the crystal growth rate and direction are given by the topology of the separated amorphous phases surrounding the crystal. In such systems, the increase of the overall crystallinity is highly irregular, depending on whether a highly concentrated amorphous zone is available around the crystal at each time. Furthermore, the effect of the interfacial growth rate on the final structure is remarkable. With a faster crystallization, the single crystals are more structured and form percolating pathways for each material with smaller lateral dimensions. Moreover, the higher the growth rate, the more amorphous areas can be found in the final structure, which can be considered as crystallinity defects. These findings might be crucial for the efficiency of solar cells, whereby the dimensions of the phases should remain in the order of some tenths of nanometers to ensure exciton dissociation, and for which pathways to the respective electrodes should be available for both electrons and holes. In a future work, this model will be improved to take into account further nucleation of the crystals as well as crystal growth anisotropy and partial crystallinity which are typically encountered for polymeric donor materials. Complementary to this work where we dealt with simple model systems, it will be used to investigate realistic OPV donor/acceptor polymer/small molecule systems. Thereby, the challenge is to obtain reasonable input parameters from experimental measurements. In particular, the phase diagrams of the systems should be identified [[15]] [ [34]] [ [35]] [ [37]] , focusing on the question whether liquid-liquid spinodal decomposition is expected.
Nevertheless, this work also shows that the knowledge of the kinetic parameters (crystal growth rate, diffusion coefficients but also nucleation rate) is very important to quantitatively investigate the morphology evolution.