Crystalline morphology formation in phase-field simulations of binary mixtures

Abstract

Understanding the morphology formation process of solution-cast photoactive layers (PALs) is crucial to derive design rules for optimized and reliable third-generation solar cell fabrication. For this purpose, a Phase-Field (PF) computational framework dedicated to the simulation of PAL processing has recently been developed. In this study focused on non-evaporating, crystallizing binary mixtures, distinct crystalline morphology formation pathways are characterized by a systematic exploration of the model's parameter space. It is identified how, depending on material properties, regular, dilution-enhanced, diffusion-limited and demixing-assisted crystallization can take place, and which associated structures then arise. A comprehensive description of the thermodynamic and kinetic mechanisms that respectively drive these separate crystallization modes is provided. Finally, comparisons with experimental results reported in the literature highlight the promising potential of PF simulations to support the determination of process–structure relationships for improved PAL production.

---

1 Introduction

With persistently increasing power conversion efficiencies (PCEs) and a broadening field of applications, solution-processed third-generation solar cells are expected to play a significant part in the transition towards more sustainable and environment-friendly energy production infrastructures. Especially photovoltaic devices manufactured from perovskite and/or organic semiconductors present promising perspectives for both low-cost and energy-efficient fabrication on an industrial scale.^1,2 Indeed, already well established solvent-based thin film deposition techniques, such as slot-die coating, are compatible with these materials and can be integrated in a roll-to-roll process.^3,4 Nonetheless, despite attractive processing features, both technologies still face numerous challenges on the path towards reliable lab-to-factory upscaling and widespread commercialization.^5,6 Although the nature of the confronted issues is specific to the selected material properties, the nanostructure of the photoactive layer (PAL) is always determinant for the entire device characteristics, and, most importantly, for PCE. Since the morphology of the PAL arises during film fabrication, understanding how it is controlled by process parameters is of paramount importance to produce performant and robust solar cells.

At present, advances regarding optimal PAL processing are still mainly achieved by following a time- and resource-costly trial-and-error approach and many process–structure relationships remain yet to be unraveled. Specialized numerical simulations are therefore of interest to gain theoretical insights into the complex interaction of physical phenomena that shapes the PAL at nanometric scale.⁷ For instance, organic photovoltaic (OPV) material systems are likely to experience a miscibility gap while undergoing temperature and composition changes during processing.^8,9 Knowing under which conditions amorphous–amorphous phase separation (AAPS) interferes with other structure formation mechanisms is thereby of crucial interest to pinpoint the best way for device fabrication. For this purpose, multiphysics phase-field (PF) modelling frameworks were employed by Wodo et al.,¹⁰ Kouijzer et al.,¹¹ and Schaefer et al.,¹² to examine the coupling of AAPS and evaporation in ternary solvent-donor–acceptor blends for OPV. Furthermore, single component crystallization in OPV films produced via meniscus-guided coating was also addressed by Michels and co-workers.¹³ Finally, a PF framework combining all these phenomenologies was recently designed to compute the time-resolved evolution of solution-cast PAL architectures over the whole fabrication process.¹⁴ In previous reports, most relevant mechanisms involved in the manufacturing of bulkheterojunction-type PALs (BHJ), namely spinodal decomposition (SD),¹⁵ its interplay with crystal growth,¹⁶ and solvent evaporation^17,18 were successfully implemented and analysed in-depth. A complete investigation of the film drying, which encompasses all aforementioned physics along with certain aspects of fluid dynamics, was also carried out for the extensively surveyed P3HT-PCBM OPV model system.¹⁹

A comprehensive study of the progress of crystallization depending on available morphology formation pathways, which are recognized to be strongly impacted by competing processes such as AAPS,²⁰ remains to be realized. Therefore, it is now sought to obtain a general overview of the crystallization based on the physics that are captured by this PF framework. More precisely, the objective is to understand how material and processing parameters influence the crystallization kinetics and the resulting morphologies. The scope is restricted to non-evaporating incompressible binary mixtures with only one crystallizing species, as this configuration readily provides numerous complex cases for systematical analysis. Note that hydrodynamic effects are also neglected. Nevertheless, the application range extends beyond the material classes used in the field of solution-processed photovoltaics. In this broader context, the PF method is a well-established continuum mechanics simulation approach that allows separately for modelling of crystallization²¹ and amorphous–amorphous phase separation.²² Different couplings between both phenomena have been derived to investigate crystallization processes in specific systems such as nucleation assisted by phase separation fluctuation in PEH–PEB polymer–polymer blends,²³ metal alloy solidification from regular and phase separating solutions²¹ or microstructure evolution in controlled drug-release coatings.^24–26 A development, which the current code relies on, has been proposed by Matkar and Kyu²⁷ to account simultaneously for crystallization and AAPS and permit supplementary crystalline–amorphous and crystalline–crystalline demixing. In this work, two main tasks are associated with the overall aim: first, a parametric study is conducted to elucidate the sensitivity of the crystallization dynamics against thermodynamic and kinetic parameter variations and reveal possible pathways for crystalline morphology formation. Second, the identified distinct crystallization modes are further categorized with respect to the underlying physical properties that activate them.

The remainder of this paper is organized as follows: Section 2 focuses on the description of the PF model while Section 3 supports the theoretical analysis of the simulated crystallization behavior with discussions about phase equilibria and potential nucleation pathways through the free energy landscape. Section 4 then reports both the morphology and the crystallization kinetics observed for the initial OPV reference simulation setup in the regular crystallization mode. Additionally, results from the sensitivity study are presented here. In Section 5, diffusion-limited, demixing-assisted and dilution-enhanced crystalline structure formation scenarios are successively characterized and noteworthy arising morphological features are detailed. A short literature survey focusing on OPV systems is undertaken in Section 6 to validate simulation outcomes against already published experimental results. Finally, the conclusions of this work are presented in Section 7.

2 Description of the phase-field model

2.1 The free energy functional

The PF approach makes use of field variables to represent the local phase properties of the investigated material system. Since the PF variables vary continuously across the simulation domain, interfaces between different phases that arise within the system are inherently diffuse. Here, with the aim to simulate crystallization in incompressible binary blends where only one component can crystallize, two different parameters denoted by ϕ and ψ are employed. On the one hand, ϕ is a compositional PF variable which monitors the volume fraction of the crystalline component. Even though the scope of this study is broader than solvent–solute mixtures, the crystalline species is abusively referred to as “solute” in what follows for simplicity. In contrast, the amorphous one is labeled as “solvent”. On the other hand, ψ is a structural order parameter which tracks the extent of the crystallization transformation. It also ranges from 0 (fully amorphous state) to 1 (fully ordered state, i.e. completed crystallization).

Depending on the values of the field variables ϕ and ψ, successive evaluations of the system's free energy dictate how it evolves in time. The PF model used in this article is a particular case of a more general previously established framework¹⁴ to focus on the coupled physical phenomena that are relevant for the present problem, namely crystallization and phase separation. The former can be defined as the transition from a micro-structurally disordered phase, termed amorphous phase, to an energetically more favorable phase with compact ordered micro-structure: the crystal. For a crystal to form, an initial localized energy fluctuation is always required to enable the system to build a stable nucleus that overcomes the crystallization energy barrier and start free crystal growth.²¹ For this reason, crystallization is categorized as a nucleation and growth (NG) process. At late crystallization stages, free growth can be hindered by crystal impingement, lack of material, inadequate diffusion properties or other boundary conditions.

The second process, phase separation, is the demixing of a blend of at least two species which mixed state is thermodynamically unstable (or metastable). As a result, distinct, well-defined domains, which are respectively richer in each of the separating components, are formed within the system. Differently from crystallization, phase separation can either be triggered locally through NG when the mixture is metastable, or happen spontaneously and simultaneously in the whole blend when it is unstable. In the latter case, the demixing is referred to as spinodal decomposition (SD). In addition, it can be mentioned that species involved in the dissociation process might be in the crystalline or in the amorphous state. In this report, amorphous–amorphous phase separation is abbreviated with AAPS.

In order to describe the thermodynamics of both phase change phenomena, the Gibbs free energy G of the system is written as the sum of two corresponding contributions,

G = G^(ac) + G^(mix), (1)

where G^(ac) represents the free energy change due to transitions from the amorphous to the crystalline state, while G^(mix) captures interactions between both mixed components that determine whether the system is prone to phase separation or not. Following the treatment by Cahn and Hilliard²⁸ or Ginzburg and Landau,²⁹ the free energy of the system is typically expressed as a functional in the context of PF modelling. The next subsections detail further how both terms of eqn (1) are formulated in the employed framework.

2.1.1 Free energy of crystallization. The free energy change due to crystallization (G^(ac) in eqn (1)) is rewritten as a functional that includes three main contributions:whereby, for the present PF simulations, V indicates the volume of the discrete mesh elements which the simulation box consists of. The first term in eqn (2) accounts for the bulk free energy of crystallization phase change and has to take the form of a double-well potential in order to mimic the behavior of NG processes. Usually, polynomial interpolation functions are exploited to provide this property.²¹ Here, a combination of a quartic function q(ψ) = ψ²(1 − ψ)², which creates the desired double-well shape, and a cubic interpolation p(ψ) = ψ²(3 − 2ψ), that regulates the latent heat release with increasing order parameter ψ, is employed. The polynomial q(ψ) is proportional to a coefficient W (in units J kg⁻¹) which is related to the crystal's surface tension and controls the height of the energy barrier to be overcome before NG. Likewise, p(ψ) is weighted by a factor L(T/T_m − 1), which stems from Turnbull's linear approximation of the free energy decrease for isotropically crystallizing species at low undercoolings.²¹ Analogously for lamellar polymer crystallization, a similar factor emerges in Lauritzen and Hoffmann's treatment.³⁰L (in J kg⁻¹ as well) represents the latent heat of fusion (assumed here to be constant), while T and T_m respectively denote the system's temperature and the melting temperature of the crystalline component. The total volume specific bulk free energy change due to crystallization finally readswith ρ the density of the species subject to crystallization.

Furthermore, the second term in eqn (2) describes the energy rise pertaining to transitions in molecular ordering at crystal boundaries:

Therefore, it scales with the squared gradient of the order parameter |∇ψ|² as well as a coefficient ε², which, similarly to W, is linked to the crystal surface tension.³¹

The last term in eqn (2) introduces a supplementary PF variable θ that attributes an orientation angle to nucleated crystals and writes

An orientation is assigned to crystallizing regions of the domain once an empirically determined value is surpassed locally by the order parameter ψ. If this threshold is attained in the vicinity of a location where θ is already defined, the same orientation is given, so that the θ-field expands according to crystal growth. Otherwise, a new θ-value is attributed (new nucleus). Conversely, if the order parameter drops below the threshold, the orientation is locally suppressed (more information about this can be found in a previous publication¹⁴). For the present study, the used criterion is ψ > 0.5. This limit is set to be higher than the order parameter for which the direct nucleation energy barrier is reached (see Section 3.2 and SI-C, ESI†). In this way, an orientation is only given when the bulk free energy for crystallization is favorable for stable crystal build-up.

The parameter θ is required to handle crystal impingement due to orientation mismatch and prevent interpenetration. An energy coefficient α scales the related energy contribution (eqn (5)). δ(∇θ) essentially works as a step function that becomes 1 when a gradient in orientation is (a) defined and (b) different from 0. Thus, the orientation mismatch energy only becomes active when two distinctly oriented crystals are on the verge of impingement. The interpolation polynomial p(ψ) then regulates the strength of G^(ori,ac)_V according to the advancement of crystallization at the impingement location. Note that θ can be used to implement anistropic crystal growth,³² which is often observed for polymers and organic small molecules. Here, however, the crystals are maintained isotropic since anisotropy is not expected to significantly impact the investigated trends in crystallization kinetics.

2.1.2 Free energy of mixing. Analogously to the crystallization term, the free energy change due to mixing properties can be split into several contributions:

Again, G^(bulk,mix)_V and G^(grad,mix)_V designate bulk and interfacial free energy densities, whereas G^(num,mix)_V is a supplemented term designed to ensure numerical stability when volume fractions reach values close to 0 or 1.^14,24 This contribution reads as

with both numerical parameters β and γ empirically set to 10⁻⁵ J m⁻³ and 1, respectively, so as to not significantly impact the physics of the simulations while still providing the sought regularization in the direct vicinity of the boundaries of the ϕ-interval.

According to the Flory–Huggins theory of mixing,³³ which allows for modelling of small molecule solutions up to polymer mixtures, and its generalization by Matkar and Kyu²⁷ to include crystalline components, the bulk contribution G^(bulk,mix)_V can further be written as

Here, R is the ideal gas constant and v₀ the molar volume of the Flory–Huggins lattice elements, which can be as large as atoms, monomers, polymer segments or small molecules, depending on the selected material system. N₁ and N₂ then express the respective species sizes in terms of these elements. Throughout this report, subscript 1 refers to the crystallizing component and 2 to the amorphous one. Since the amorphous species always possesses the smallest molar volume in the present investigations, v₀ is indexed on it and N₂ is equal to 1.

χ _aa denotes the so-called Flory–Huggins interaction parameter which accounts for pairwise interactions between neighbouring Flory–Huggins lattice elements in the amorphous state. The nature of these interactions, which determine the miscibility of the amorphous mixture, may be enthalpic or entropic. In general, a linear relationship in 1/T is assumed for the temperature-dependence of χ_aa:³⁴

The offset A is associated with entropic contributions while the slope B sets the enthalpic interaction trend in 1/T. In experiments, different methods are commonly employed to evaluate χ_aa.^11,35–37 Nevertheless, adjustment of the classical Flory–Huggins model was suggested by Matkar and Kyu²⁷ to account for separate interaction contributions that take place between components in the amorphous and the ordered state. To this end, an auxiliary interaction parameter χ_ca that satisfies the following relation was introduced:

In this expression, the constant C is a proportionality factor. It can be mentioned that this is in qualitative agreement with the Flory–Huggins theory modification by Hu³⁸ that also results in a corrected Flory–Huggins interaction parameter with a supplementary term proportional to a crystallization energy. Additionally, the multiplication of χ_ca with ψ² in eqn (8) couples ordering transitions to composition changes, so that crystallization of one species is likely to be accompanied by a demixing process. Besides temperature, χ_aa and χ_ca are also expected to change with blend composition.³³ For simplicity, they are however assumed independent of ϕ in this work.

Finally, like the order parameter ψ, gradients in volume fraction ϕ at phase interfaces give rise to an energy contribution related to surface tension. The corresponding term, G^(grad,mix)_V in eqn (6), scales with a coefficient κ which is comparable to ε² in G^(grad,ac)_V (eqn (4)):

2.2 Kinetic transport equations

2.2.1 Stochastic Allen–Cahn equation. In order to track the time-evolution of non-conserved order parameter fields, the Allen–Cahn equation is usually employed.³⁹ For the purpose of this work it is augmented with a Langevin force term to reproduce the effect of thermal fluctuations on the system. Under this form, it is known as stochastic Allen–Cahn equation and reads asHere, δG/δψ symbolizes the functional derivative of the free energy with respect to ψ, which, once written in terms of the corresponding partial derivatives of the volume specific free energy, yieldsM(ϕ) denotes a kinetic parameter that determines the local conversion rate of amorphous material into crystalline one. Therefore, it is determinant for crystal nucleation properties and interface mobilities during growth. As it is further detailed in an upcoming section (Section 5.4), M is expected to be composition-dependent²⁴ in the general case. Finally, the Langevin term ξ is a Gaussian noise with a standard deviation of designed to preserve the fluctuation–dissipation theorem.⁴⁰ It is responsible for thermally activated nucleation, as well as crystal coarsening.

In principle, an Allen–Cahn-type equation similar to eqn (12) can be used to describe the dynamics of the orientation θ. The orientation of crystals may change over time, which can generate grain melting or crystallographic branching in anisotropically crystallizing materials that form dendrites and/or spherulites.⁴¹ In this study, the crystals are isotropic and the orientation field is only used to differentiate between crystals and handle their impingement. For this reason, no additional differential equation needs to be solved for, the evolution of the θ-field (i.e. apparition and suppression of an orientation) being handled by the order parameter threshold, as already reported in Section 2.1.

2.2.2 Cahn–Hilliard conservation equation. Conserved quantities such as the volume fraction ϕ satisfy a Cahn–Hilliard-type continuity equation⁴² which equates their partial derivative with respect to time to the divergence of their corresponding flux J. Considering assumptions from non-equilibrium thermodynamics theory for multicomponent systems,⁴³J can be substituted as follows:

The quantity μ₁ − μ₂ is a chemical potential density which is sometimes referred to as exchange chemical potential^9,11 and which derives from the classical diffusion potential.^18,44 Its gradient is the thermodynamic driving force for mass diffusion. According to Cahn,⁴⁵ it can be computed as the functional derivative of the free energy with respect to ϕ, i.e.

Furthermore, Λ(ϕ) designates the Onsager mobility coefficient,⁴⁶ which, in the more general multicomponent case, takes the form of a positive semi-definite matrix. Diffusion of species within the system is regulated by this kinetic parameter. To calculate Λ depending on ϕ, the PF code implements the fast⁴⁷ and slow⁴⁸ mode theories. Both have been confirmed to lack generality so that attempts to unify them have been published.⁴⁹ However, the fast mode is selected here, as a strong impact of the chosen theory on the qualitative trends described in this article is not expected. Λ(ϕ) accordingly writes as

Λ(ϕ) = (1 − ϕ)²ϕN₁D^(self)₁(ϕ) + ϕ²(1 − ϕ)N₂D^(self)₂(ϕ), (16)

where D^(self)₁(ϕ) and D^(self)₂(ϕ) stand for the self-diffusion coefficients of species 1 and 2, respectively. To interpolate their values in the volume fraction range from measurements in pure materials (D^(self)_1|ϕ→0, D^(self)_2|ϕ→0, D^(self)_1|ϕ→1, D^(self)_2|ϕ→1), a logarithmic mean is utilized for simplicity (see SI-A, ESI†).

Note that other dependencies have been suggested, as estimating the composition-dependence of self-diffusion coefficients is still an open research topic.⁵⁰ Within crystals where diffusion is expected to be slower by several order of magnitudes than in the amorphous state, an additional penalty function decreases the coefficients D^(self)₁(ϕ) and D^(self)₂(ϕ) depending on the value of the order parameter ψ. More details about this can be found in previous publications.^14,16,19 For this work, the parameters of this penalty function are taken identical to those used in a precedent study by Ronsin and Harting.¹⁹

Similarly to the stochastic Allen–Cahn transport equation, fluctuations can be added in the Cahn–Hilliard formula (eqn (14)), which is then referred to as Cahn–Hilliard–Cook equation.⁵¹ These fluctuations are mainly of use to trigger AAPS. It was noticed during preliminary experiments that fluctuations in the stochastic Allen–Cahn equation (eqn (12)) already serve this purpose due to the coupling between crystallization and phase separation free energies (eqn (8)). The noise contribution on the Cahn–Hilliard equation is therefore omitted to ensure better numeric robustness, especially regarding volume fraction conservation.

3 Theoretical analysis

3.1 Phase equilibria

In accordance with Gibbs' phase rule, binary systems under constant pressure and temperature conditions can present at most two distinct phases in the equilibrium state. Since the PF model relies on a free energy formulation, information about the phase equilibrium compositions can already be obtained before carrying out numerical simulations of the system's evolution in time. Here, the approach described by Horst^52,53 is employed in order to produce ϕ − T phase diagrams. Note that, although solely this sort of graphs is featured in this article, this routine can as well be used to trace figures with different ordinate axes, as for instance ϕ − χ_aa diagrams, which are also often examined in the literature.^{9,15,36,54,55}

The typical form of the phase diagram corresponding to systems investigated in the present work can be visualized in Fig. 1. The increasing divergence of the liquidus and the solidus from the pure solute melting temperature towards lower T is notably due to the UCST-type (upper critical solution temperature) behavior induced by χ_aa (with positive coefficient B in eqn (9)). Even though multitudes of diagram shapes, including LCST-types (lower critical solution temperature), can be achieved with this free energy,⁵⁶ this sort of diagram is considered as a well-suited starting point for the parametric study since it is frequently reported in the OPV literature.^9,57–59

Fig. 1 Phase diagram for the PCBM-oDCB system as extrapolated with the current free energy model from experimental solubility measurements by Schmidt-Hansberg et al.⁶⁰ (χ_aa = 0.34 + 128.14/T and χ_ca = 0.1648 v₀N₁ρL/RT). Other thermodynamic parameters used for the computation are listed in Table 1.

The reference parameter set that corresponds to the phase diagram depicted in Fig. 1 represents a mixture of (6,6)-phenyl C61 butyric acid methyl ester (PCBM) dissolved in ortho-dichlorobenzene (oDCB). This OPV small molecule–solvent blend is selected as a practical initial system due to the availability of experimental data from which inputs for the PF model can be derived.¹⁹ In the OPV community, PCBM is well-known as the electron acceptor of the extensively studied P3HT-PCBM polymer-small molecule bulkheterojunction (BHJ).^59,61,62 More recently, PCBM has also been mixed with non-fullerene acceptors (NFAs) in ternary and quaternary OSCs to improve device performance.⁶³Ortho-dichlorobenzene (oDCB) is a high-boiling point solvent commonly used for its affinity with a wide variety of OPV materials.⁸

Thermodynamic parameters employed for the calculation of PCBM-oDCB phase equilibria are summarized in Table 1 and have been discussed in detail in previous numerical investigations of this system.¹⁹ In this work, the energy barrier coefficient W is slightly modified to grant the necessary double-well shape of the crystallization free energy over a sufficiently large parameter range for the sensitivity study (see limiting criterion on W in SI-C, ESI†). Moreover the interaction parameters χ_aa and χ_ca are adapted to fit solubility limit measurements published by Schmidt-Hansberg et al.⁶⁰ A qualitative summary of the impact of thermodynamic parameter variations on liquidus and solidus compositions is given in the ESI† (SI-B). It can be seen that an increase in both χ_aa and χ_ca shifts the position of the liquidus in opposite directions along the ϕ-axis. Therefore, several [χ_aa, χ_ca] pairs could be used to extrapolate virtually valid diagrams from the experimental data points. The availability of the exact PCBM-oDCB phase diagram is not a mandatory requirement for the present endeavor. The initial parameter set is thus selected so that enough margin is provided for the parameter variations of the sensitivity study and to avoid numerical convergence issues.^14,31

Table 1 List of thermodynamic and kinetic parameters used for the initial reference simulation of the PCBM-oDCB binary blend. The scaling parameter for the self-diffusion coefficients (Λ₀) is set equal to 1 in all simulations for which no other value is explicitly specified (see Sections 4.2, 5.1 and 5.3)

Parameter	Value	Units
ρ	1600	kg m^−3
v 0	1.131 × 10^−4	m^3 mol^−1
N 1	5.0298	—
N 2	1	—
T	333	K
T m	558	K
L	20 000	J kg^−1
W	40322.5806	J kg^−1
χ aa	0.7248	—
χ ca	1.0836	—
ε ^2	1 × 10^−10	J m^−1
κ	2 × 10^−10	J m^−1
α	8.1621 × 10^7	J m^−3
M 0	0.1	s^−1
Λ 0	1	—
	5Λ0 × 10^−10	m^2 s^−1
	Λ 0 × 10^−13	m^2 s^−1
	2Λ0 × 10^−9	m^2 s^−1
	Λ 0 × 10^−12	m^2 s^−1

---

3.2 Nucleation pathways through the free energy landscape

The analysis of phase diagrams allows to predict the equilibrium states towards which thermodynamic systems eventually evolve. However, these diagrams contain no information about intermediate states which are explored on the way to the final equilibrium. This renders it difficult to ascertain under what conditions nuclei effectively appear and grow. In order to gain further understanding of the privileged routes for phase change at a given process temperature, the topology of the corresponding free energy surface can be examined. Fig. 2a shows an example of such a three-dimensional free energy landscape computed for the PCBM-oDCB small molecule–solvent blend at 60 °C. In the following, crystalline germ formation is considered from an originally amorphous blend (i.e. ψ = 0) with homogeneous composition ϕ₀. The trigger for crystal nucleation is the stochastic term ξ in the Allen–Cahn equation (eqn (13)), which is responsible for the emergence of regions with a non-zero order parameter ψ in the mixture.

Fig. 2 (a) Three-dimensional free energy landscape computed for the PCBM-oDCB system at T = 333 K. (b) Corresponding 2D contour plot with highlighted demixing and crystal growth regions.

The first possible nucleation pathway through the landscape allows for stable crystallites to be produced in a single step. More precisely, if, over an area wider than the critical nucleus size, fluctuations on ψ arise large enough to overcome the energy barrier, then crystal ordering (i.e. an increase in ψ), growth and purification in terms of solute content (i.e. an increase in ϕ) can spontaneously proceed. The second pathway is a two-step process which involves a demixing stage after the initial order parameter fluctuation and before the stable crystal build-up. To understand when each nucleation route is followed, the contour plot of the free energy landscape (Fig. 2b) can be used. In Fig. 2b, the red and black lines respectively stand for the maximum (i.e. the energy barrier) and minimum of the free energy at fixed solute volume fractions (see SI-C for the corresponding analytical formulae, ESI†). The final, global equilibrium points of the thermodynamic system, which correspond to the liquidus and solidus in the phase diagram, are reproduced in purple.

The one-step pathway is only available if the crystallization occurs in a mixture with a blend ratio ϕ₀ for which an energy barrier exists (in the example of Fig. 2b this corresponds to ϕ₀ > 0.2). There, if, over a domain that is larger than the size of the critical nucleus, simultaneous local fluctuations surpass the energy barrier (see Section 2), a crystal germ can form directly. Moreover, if, beyond the barrier, the minimum associated with the overall composition ϕ₀ is a global minimum (solid black line in Fig. 2b), the germ is stable and can only evolve towards the solidus point. This occurs by further crystal ordering and purification (i.e. increases in ψ and ϕ, respectively). In case the minimum after the barrier is local (dashed black line in Fig. 2b), the resulting germ relies on composition changes due to crystal purification to stabilize. Therefore, if crystal purification is hindered for example due to strong diffusion limitations (see Section 5.1), one-step crystallization can be kinetically impeded below a threshold volume fraction (here ϕ₀ ≤ 0.4).

In comparison, two-step nucleation can take place when the initial fluctuation locally drives the system into a region of the landscape where the composition ϕ₀ is energetically unstable or metastable (represented by the solid and dashed golden lines in Fig. 2b, see SI-C for analytical derivation, ESI†). The demixing is undergone in a phase dissociation process similar to AAPS with spinodal- and binodal-like behavior. However, unlike common AAPS, the unstable (or metastable) state is in general not reached in the whole system due to the local nature of the order parameter fluctuations. Thus, the demixing is always inhomogeneous and tends to proceed by nucleation and growth. It can be remarked that the regular AAPS in immiscible amorphous blends is in fact a particular case of this more general configuration where the strength of the fluctuations required to enter the demixing region is minimal. A free energy landscape for a system subject to a regular homogeneous AAPS can be visualized in Fig. 8b of Section 5.2. While the pseudo-AAPS develops towards the compositions predicted by the binodal-like curves, crystallization can start as soon as the region delimited by the energy barrier and the minimum (red and black lines in Fig. 2b) is entered. The exact volume fraction at which nucleation starts then depends on the balance between the kinetics related to the Cahn–Hilliard (eqn (14)) and Allen–Cahn (eqn (12)) equations.

For the particular PCBM-oDCB reference system, for which the landscape is shown in Fig. 2a and b, it can be noted that the energy barrier is located at relatively low ψ values and spans over the majority of the ϕ-range (ϕ ∈ [0.2, 1.0]). In contrast, the AAPS-like region requires high fluctuations on ψ to be attained. The two-step pathway is therefore much less probable than the one-step counterpart in this case. Below ϕ ≃ 0.2, where solely two-step processes can lead to crystallization, the appearance of nuclei is additionally strongly delayed. Indeed, at these lower overall solute concentrations, the solute majority phase which emerges from the demixing, and which reaches the landscape region that is beyond the energy barrier, is likely to be smaller than the critical nucleus size, so that, in the end, the system tends to collapse back to the original state of the mixture.

On top of these considerations, additional comments can be made: first, multiple different free energy landscape topologies can be obtained depending on the thermodynamic parameter set. Accordingly, the AAPS-like region in Fig. 2b can overlap in various manners with the energy barrier curve and modify the volume fraction ranges dominated by each pathway. For instance, two-step nucleation is favoured over a significant composition range for the immiscible blends in Sections 5.2 and 5.3. In any case, the amplitude of the fluctuations has to be sufficient to overcome the energy barrier threshold at some point of the process. In particular, this generally limits one-step nucleation from lower solute volume fractions since the associated curve increases with lower ϕ (see Fig. 2b). Second, as already mentioned, nucleation via a two-step process can be critically hindered at the lowest concentrations. As a matter of fact, since the proportions of the intermediate phases respect a lever rule, it can become very difficult to start from low ϕ₀ and generate domains which are simultaneously larger than the critical nucleus size and have an adequately high solute content. Third, kinetic effects also play a crucial role. In general, the exact path through the free energy landscape is determined by the respective timescales on which the transport phenomena occur. In extreme cases where the evolution of the volume fraction is very slow as compared to the order parameter ψ, two-step nucleation can be significantly impaired (see Section 5.1). Then, the system regresses rapidly to the amorphous state (ψ = 0) before AAPS-like demixing produces phases with favorable compositions for nucleation of stable crystals.

4 General crystallization behavior in binary mixtures

4.1 Regular crystalline morphology formation and reference transformation kinetics

Following these preliminary discussions, simulations of crystalline morphology formation in binary mixtures are now addressed. For this purpose, all computations carried out in the context of this work rely on a 2D square box of 512 nm² with a grid spacing of Δx = Δy = 1 nm as simulation domain. Periodic boundary conditions are applied at each edge to model crystallization within the bulk of the blend. In order to optimize time-stepping, the time scheme makes use of an adaptive procedure that allows Δt to reach the largest possible value which still ensures numerical convergence.¹⁴

To begin the exploration of the different possible scenarios, the crystallization regime herein referred to as “regular” is considered, as it produces the most similar structures in comparison with crystallization from a pure material melt. In this configuration, the binary blend is miscible in the amorphous state and solute material is always available at the crystal growth front due to its relatively rapid diffusion towards interfaces. This means that separation of both chemical species only occurs once a crystal germ has formed and that the growth of a crystal is controlled by the intrinsic mobility of its interface, which is determined by M(ϕ). In addition, the inclusion rate of molecules into a crystal is assumed to be independent of the local mixture composition, i.e. M(ϕ) = M₀. This constant mobility parameter M₀ should be calibrated with experimental crystallization kinetics data when specific material systems are to be replicated. Here, it is initially set to 0.1 s⁻¹, as this value (in combination with the other initial parameters) is suitable to realize the regular crystallization mode. The self-diffusion coefficients necessitated for calculation of Λ(ϕ) (see eqn (16) and (17)) are taken over from earlier studies on the PCBM-oDCB solution.¹⁹ In addition, the value of the impingement energy coefficient α is selected so that crystal interpenetration is prevented. The full set of thermodynamic, kinetic and numerical parameters used for the computations is summarized in Table 1.

Fig. 3 compares the progress of crystallization for three different initial blend ratios (ϕ₀ = 0.4, ϕ₀ = 0.6, ϕ₀ = 0.8). First, it can be observed that, for every simulated setup, three stages of crystallization, namely nucleation, free growth, and grain coarsening – i.e. the expansion of bigger crystals at the expense of smaller ones – are present. Apart from the final amount of crystallized material, which is straightforwardly related to the initial overall content of solute ϕ₀, the simulations mainly differ in the extent in which these different mechanisms shape the crystallization process. More specifically, the blend ratio affects the balance between nucleation and growth. In this way, it can be remarked that morphology formation from lower ϕ₀ tends to be more growth-dominated than from higher blend ratios – which, in turn, is nucleation-dominated. As can be visualized on the free energy surface contour plot (Fig. 2b), the value of ψ which corresponds to the location of the energy barrier increases with lower ϕ. For this reason, less fluctuations are able to overcome the barrier. The associated critical nucleus size augments as well. As a consequence, stable nuclei are less likely to appear, leaving more solute material available for incorporation in already formed germs. In this configuration, the growth stage is favored compared to nucleation and only a limited number of relatively large crystals emerges. Conversely, with higher solute content, the energy barrier location drops to lower ψ. Accordingly, nucleation is facilitated up to the point where growth is hindered since most of the crystalline material is directly consumed for nuclei production. Crystals are then comparatively smaller (at least, before coarsening becomes significant).

Fig. 3 Progress of crystallization for the PCBM-oDCB binary system at T = 333 K and for three different blend ratios ((a) ϕ₀ = 0.4, (b) ϕ₀ = 0.6, (c) ϕ₀ = 0.8). The volume fraction field ϕ of the crystallizing component is shown with a color code ranging from deep red (high ϕ) to deep blue (low ϕ). Thermodynamic and kinetic parameters used for the calculation are specified in Table 1.

Grain coarsening occurs predominantly at the end of the crystallization process, when most of the available solute material has effectively crystallized. It is driven by the minimization of surface tension contributions to the free energy of the system, which are most important at grain interfaces. This leads to a reduction of boundaries in between crystals. Remark that this coarsening mechanism is mathematically expressed by the gradient terms in eqn (13). Another coarsening phenomenon that acts on phases of distinct compositions is captured by eqn (14) and will be more visible in an upcoming section (Section 5.2). From a general perspective, this shows that the proportions in which the species are initially mixed can lead to distinct formation regimes and eventually produce different structures with variations in crystal size distributions and spatial arrangements.

The crystallization kinetics can now further be characterized by recording the time-dependent evolution of the extent of transformed material relatively to the total amount of solute. This quantity is labeled as “crystallinity” in what follows and is compared for several initial blend ratios in Fig. 4a. It can first be noticed that all curves present a sigmoidal shape, which is frequently observed for crystallization processes and usually modelled with the Johnson–Mehl–Avrami–Kolmogorov (JMAK) theory^64,65 and its extensions.⁶⁶

Fig. 4 (a) Transformation kinetics in the regular crystallization regime at different PCBM-oDCB blend ratios. (b) Corresponding crystallization half-time τ_50% as a function of blend composition.

Moreover, it can be seen that the higher the global volume fraction of solute, the faster is the overall transformation. This is expected since the thermodynamic driving force for crystallization is proportional to ϕL(T − T_m)/T_m (see eqn (3)). It is worth pointing out that, in this configuration, the coefficient for the local crystallization rate is kept constant (M(ϕ) = M₀). For systems where diffusion properties vary significantly with the blend ratio, a related composition-dependence of M(ϕ) is anticipated to impact substantially the crystallization kinetics with different ϕ₀. Further details about this are given in Section 5.4.

Finally, recalling the discussion about nucleation pathways through the free energy landscape (Section 3.2), it could be verified that no crystallization occurs within reasonable computation times at the lowest ϕ₀. The reported volume fraction range thus starts at 0.2 in all forthcoming figures. Besides, even for higher mixing ratios that allow for direct one-step processes, the onset of crystallization can be considerably delayed because of the importance of the surface tension contributions (i.e. ε², and W) that impede nuclei formation. Thus, no crystallization is witnessed below ϕ₀ = 0.4 within the timespan simulated in Fig. 4.

4.2 Thermodynamic and kinetic sensitivity study

Having analysed the crystallization behavior over the whole concentration range for a fixed binary system (Table 1), the sensitivity of the process kinetics against thermodynamic and kinetic parameter variations is now investigated. To quantify this, the half-time of crystallization τ_50% is defined as the time when the extent of crystallized material reaches 50% of the total amount of available crystalline component. τ_50% can be extracted from transformation plots such as Fig. 4a. The overall crystallization rate can then be approximated by taking its inverse (i.e. 1/τ_50%). In order to carry out this study, all parameters from Table 1 are kept constant, except the indicated variable of interest, which is examined within the interval reported in Fig. 5.

Fig. 5 Sensitivity of the crystallization half-time τ_50% to (a) and (b) surface tension, (c)–(f) thermodynamic and (g) and (h) kinetic model parameter variations for different initial blend ratios ϕ₀. For each separate benchmark, all parameters are maintained constant at the values reported in Table 1, except the variable specified on the abscissa.

At first, the focus is restricted to surface tension parameters, namely ε² (eqn (4)) and κ (eqn (11)). As can be seen in Fig. 5a and b, ε² has globally a more severe influence on the crystallization kinetics as compared to κ. In addition, the general sensitivity is more important at low volume fractions for both coefficients. It can be remarked that the curves for the more growth-dominated low ϕ₀ tend to fluctuate more due to the emergence of fewer nuclei, which results in poorer statistics for τ_50%. It is also noteworthy that the balance between both parameters affects which of the surface tensions arising at interfaces between a crystalline region and its amorphous surroundings, or at boundaries separating two distinct crystal grains, is the strongest. Accordingly, crystals tend to agglomerate or spread across the domain, depending on which associated energy contribution is the highest.

Next, the thermodynamic parameters L, W, χ_ca and χ_aa are varied. From Fig. 5c it can be concluded that increasing the latent heat of fusion L boosts crystallization. An increase of the energy barrier coefficient W conversely inhibits it (Fig. 5d). This is expected since the former parameter determines the thermodynamic driving force for crystallization while the latter dictates the height of the energy barrier to be overcome for nucleation. Considering the logarithmic axes of the figures, it can be seen that both have an impact on the crystallization half-time that is stronger than an exponential trend. Notice also that the dependencies displayed here for L do not involve the coupling with χ_ca (eqn (10)), that is χ_ca is kept constant. Although the spacing between the curves for different mixture compositions are not identical, the general behavior is nonetheless similar when χ_ca is assumed to vary with L according to eqn (10) (see SI-D, ESI†).

Concerning the evolution of τ_50% as a function of the amorphous–crystalline interaction parameter, it can be observed in Fig. 5f that the half-time of crystallization is delayed with ascending χ_ca. This is because the energetical cost due to enthalpic interactions within crystals is larger. It can be recognized that the effect is more pronounced for lower ϕ₀. In contrast, the crystallization rate increases with χ_aa (Fig. 5e). This is coherent in terms of energy balance since strong interactions between both species in the amorphous phase promote crystallization. Similarly to the results for χ_ca, the impact is more significant at lower solute content. Past a certain χ_aa threshold, crystallization is furthermore abruptly accelerated and becomes even faster than for some of the higher concentrated solutions. This behavior is symptomatic for an AAPS occurring ahead of the crystallization process and will be detailed in an upcoming section (Section 5.2). To predict the appearance of the amorphous miscibility gap, the following analytical formula can be used:³⁴

where represents the critical value of the amorphous–amorphous interaction parameter beyond which the system might become immiscible depending on the blend ratio ϕ₀. N₁ and N₂ also play a role in this equation. In the present context, N₂ is always equal to 1. Therefore, increasing the relative solute size N₁ in the Flory–Huggins lattice reduces miscibility as well and can eventually trigger AAPS before crystal nucleation (see SI-D, ESI†).

In the end, the influence of the kinetic parameters that account for molecular ordering and diffusion rates, respectively M₀ and Λ(ϕ) (eqn (13) and (14)), remains to be analysed. Since Λ is a function of ϕ (eqn (16)), it is scaled homogeneously by multiplication with a non-dimensional coefficient Λ₀. On the one hand, it is visible in Fig. 5g that the crystallization half-time is inversely proportional to M₀ until a certain value (M₀ ∼ 10⁴ s⁻¹) beyond which the curves slightly deviate from the trend, especially for lower solute content. On the other hand, τ_50% is almost insensitive to variations in diffusion properties above Λ₀ ∼ 10⁻⁵ (Fig. 5h). Below this point the crystallization kinetics experience a slowdown, which is also more pronounced at lower ϕ₀.

Considering that the crystallization kinetics are bound simultaneously by the maximum speed at which molecules can attach to the surface of growing crystals (determined by M₀) and the rate at which solute material can diffuse from amorphous regions to crystal interfaces (given by Λ(ϕ)), the ratio of both parameters, M₀/Λ(ϕ), establishes the governing equation that drives the crystallization process. When M₀/Λ(ϕ) is low enough, diffusion in the amorphous phase occurs on timescales smaller by several orders of magnitude in comparison to crystal growth, so that there is always sufficient material available at the growth front. Then, the crystal interface mobility M₀ is the limiting kinetic factor that dictates the crystallization rate. Accordingly, a variation in Λ₀ has no effect on the crystallization half-time in this case. Once this condition is not fulfilled anymore, i.e. at larger M₀/Λ(ϕ), crystal growth becomes diffusion-limited. Morphology formation properties specific to this crystallization regime will be examined in Section 5.1. Note that in Fig. 5g crystallization is not recorded for ϕ₀ < 0.4 even when M₀ increases. This is due to the limiting criterion for the time-step size which is inversely proportional to M₀.¹⁴ Thus, the maximum times in reach with reasonable computational effort also decrease with increasing M₀, which renders it difficult to access the timescales where the first nuclei appear. The τ_50% curves for these compositions are nonetheless expected to be parallel to the other ones (at least before the deviation from the inversely proportional trend).

Finally, it can be pointed out that the effects of all parameters related to the Cahn–Hilliard equation (eqn (14)), and thus to volume fraction changes (i.e. χ_aa, χ_ca and κ), become very limited for high solute volume fractions ϕ₀ since, in that case, the composition remains nearly homogeneous over the transformation process. Moreover, it may be mentioned that nucleation is, in general, more heavily dependent on thermodynamic parameters than crystal growth. As a result, most of the related τ_50% sensitivities presented above reflect predominantly the impact of parameters on the nucleation rate. Thus, this means that, in addition to being generally favoured at low solute content, crystal growth is also more dominant when nucleation is penalized due to high ε, W, χ_ca, or conversely low L and χ_aa values.

5 Review of noteworthy crystalline morphology formation regimes

5.1 Diffusion-limited crystallization

In this section, the morphology formation is analysed for diffusion-limited systems. In the regular crystallization case, where M(ϕ)/Λ(ϕ) is sufficiently low, diffusion within the amorphous phase is faster than crystal growth. As can be visualized in Fig. 3 (Section 4.1), this translates into a monotonously increasing volume fraction profile from the amorphous region to the crystalline one. In contrast, when molecular diffusion happens on larger timescales than crystal growth, i.e. comparatively higher M(ϕ)/Λ(ϕ), solute material is missing at the growth front since crystals consume it faster than it can diffuse from the bulk of the amorphous phase to the interface. As is shown in Fig. 6, a depletion zone consequently emerges around growing crystals. The appearance of this feature indicates the onset of the diffusion-limited regime, sometimes also referred to as diffusion-controlled. The first repercussion of diffusion limitations on the transformation kinetics is the significant slowdown of the crystallization rate noticed in Fig. 5h (Section 4.2 and SI-E, ESI†). As mentioned before (Section 4.1), systems with lower initial solute concentration are more growth dominated. Thus, and since diffusion-control affects mainly the growth stage, a decrease of the crystallization rate becomes more apparent towards low ϕ₀. Additionally, the fact that crystalline material is inherently more scarce at low overall solute contents results in more pronounced concentration deficits at the growth front, which also contributes to this kinetic hindering of crystallization.

Fig. 6 Progress of crystallization simulated for a diffusion-limited binary blend at T = 333 K. The volume fraction field ϕ is shown at different times for (a) ϕ₀ = 0.4 and (b) ϕ₀ = 0.5. Diffusion rates within the blends are slowed down by a factor Λ₀ = 10⁻⁶ in comparison to the reference simulation case (Fig. 3). Other thermodynamic and kinetic parameters are identical to those listed in Table 1. The onset of the diffusion-limited regime is characterized by the appearance of depletion zones around growing crystals. At late crystallization stages, bead chain-like structures can also be observed.

Diffusion limitations impact not only the crystallization kinetics, but also the morphology formation within the blend. Crystals tend to become round in this regime: nucleation in the vicinity of a growing germ being disabled due to the absence of sufficient crystalline material, space is available to adopt the most energetically favorable shape. It can be noticed as well that the depletion zone eventually disappears after the termination of crystal growth, when the equilibrium compositions are reached homogeneously in both amorphous and crystalline phases. Moreover, if the surface tension properties are such that boundaries between adjacent crystals are energetically favoured over interfaces from ordered to amorphous phases, crystals are prone to stick together and bead chain-like structures can ultimately be observed in the late crystallization stages, as displayed in Fig. 6.

By increasing further the ratio M(ϕ)/Λ(ϕ) the volume fraction profile progressively reacts with more delay to changes in order parameter ψ. In addition to depletion zones around growing germs, the volume fraction within crystallites becomes inhomogeneous (Fig. 7a). Indeed, it remains lower than the solidus equilibrium composition. Since a considerable amount of the second material is trapped within these impure crystals, the amount of amorphous phase is significantly lower. As a result, crystals also cover larger portions of the domain and adopt more polygonal forms. Upon coarsening, the grains become rounder but solute material inside crystals accumulates near their border rather than in the center. Although simulations were not pursued to larger timescales, it is expected that this conformation eventually breaks down and coarsens since the energetically optimal state is composed of one single crystalline droplet surrounded by the amorphous phase. The volume fractions of the phases at equilibrium should then be homogeneous, as predicted by the phase diagram (Fig. 1).

Fig. 7 Progress of crystallization simulated for a diffusion-limited binary blend at T = 333 K. The volume fraction (ϕ) and order parameter (ψ) fields are shown at different times for ϕ₀ = 0.5. Diffusion rates within the blend are slowed down by a factor (a) Λ₀ = 10⁻⁷ and (b) Λ₀ = 10⁻⁹ in comparison to the reference simulation case (Fig. 3). Other thermodynamic and kinetic parameters are identical to those listed in Table 1. In (a), the relative strength of the diffusion limitations leads to heterogeneous composition profiles within forming crystals, which consequently adopt polygonal shapes. In (b), molecular inter-diffusion is substantially slower than solute ordering. As a result, impure crystals form with a composition that remains at the initial blend ratio ϕ₀.

Finally, the diffusion-limited regime with depletion zones cannot be sustained for the highest M(ϕ)/Λ(ϕ). If the crystal interface mobility becomes too important relatively to solute diffusion, information about order parameter field fluctuations does not propagate through the volume fraction field on time, so that the composition homogeneously remains at its initial level. Two situations can be distinguished under these conditions. At lower overall content of crystalline material ϕ₀, stable one-step nucleation is thermodynamically impossible without a change in composition (see Fig. 2b in Section 3.2). Two-step processes, along with common crystal purification, are kinetically hindered as well since AAPS-like demixing follows the diffusion kinetics. The system therefore stays in its original state without crystal formation for relatively long times. As a result, it can be seen in Fig. 5g that no more crystallization is observed for low ϕ₀ (here ϕ₀ < 0.4) at the highest M₀ (or conversely at the lowest Λ₀ in Fig. 5h). In contrast, at higher concentrations, one-step nucleation is permitted (however nearly without any composition alterations due to the aforementioned kinetic restrictions). ψ values locally overcome the energy barrier and crystals form, but they remain highly impure (Fig. 7b). Equilibrium compositions and phase volume proportions can subsequently be achieved upon further grain boundary coarsening and crystal purification on considerably larger timescales.

5.2 Demixing-assisted crystallization

Crystalline structures arising in immiscible mixtures that allow for a parallel AAPS process are studied in this section. All simulations shown here base on the thermodynamic parameters used to model the PCBM-oDCB system (Table 1), by exception of χ_aa which is greater than for the investigated temperature (T = 333 K). The corresponding phase diagram and free energy landscape are reproduced in Fig. 8a and b, respectively. The unstable and metastable domains of the AAPS-like demixing region of the landscape intersect with the ϕ-axis for ϕ ∈ [0.146, 0.537] and ϕ ∈ [0.055, 0.678]. The endpoints of these intervals correspond to the classical amorphous–amorphous spinodal and binodal limits (at 333 K) which are also reported in the phase diagram (Fig. 8a), revealing a UCST amorphous miscibility gap. For blend ratios that lie within these ranges, accelerated nucleation can be triggered if the time required for AAPS¹⁵ is shorter than the time for crystallization from the homogeneous mixture at ϕ₀. In this case, amorphous demixing takes place first and generates domains with a solute content higher than ϕ₀, as expected from the binodal line. As a consequence, crystallization preferably occurs in this solute majority phase at the rate – and with the nucleation/growth balance – that is associated with this latter composition. Due to the drastic reduction in crystallization time when AAPS develops beforehand crystallization (see Fig. 5e in Section 4.2 and SI-F, ESI†), this regime is qualified as “demixing-assisted” in what follows.

Fig. 8 (a) Phase diagram of an amorphous–crystalline system which exhibits an amorphous–amorphous immiscible UCST behavior below T ≃ 440 K. The Flory–Huggins interaction parameter used for the calculation is given by χ_aa = 0.34 + 307.96/T. Other relevant thermodynamic parameters are identical to those referenced in Table 1. (b) Contour plot of the corresponding free energy landscape at T = 333 K and χ_aa = 1.2648.

In this crystallization mode, substantially different types of morphologies develop as compared to the cases investigated in the precedent sections. Fig. 9a depicts a typical formation pathway for this scenario in a non-diffusion-limited mixture (Λ₀ = 1 and M₀ = 1000 s⁻¹), i.e. the growth of crystals is controlled by the intrinsic mobility of their interfaces. In the early stages of the process, regular homogeneous AAPS (here SD) is triggered followed by coarsening. Hence, domains produced by the initial phase separation purify and tend to arrange as droplets of one component majority phase within a background matrix of the other. Volume proportions of these intermediate phases can be computed according to the lever rule applied to the amorphous binodal points of the phase diagram. In a second step, nuclei materialize in the domains of high solute content. The appearance of crystals during the coarsening partially quenches the AAPS pattern, so that it strongly determines the achieved crystalline structure. In the provided example, crystallization takes place in the background matrix, when wormlike elongated droplets of the solvent majority phase are still present. This results in the pseudo-bicontinuous conformation that can be seen at intermediate stages in Fig. 9a. Note however that multiple other geometries are possible, depending on the blend ratio, the SD or NG nature of the demixing, the quench depth for AAPS (i.e.¹⁵) and the coarsening time before nucleation. Additional cases where the blend ratio is varied can be visualized in the ESI† (SI-G). Once crystals grow, further purification of the domains in terms of solvent and solute material begins since the final equilibrium compositions are given by the liquidus and solidus in the phase diagram (blue and purple lines on Fig. 8a). A change in phase volume proportions is also predicted by the lever rule between these latter curves. As a result, it can be noticed that the crystallized solute majority domains cover less space than the amorphous background matrix they are issued from. In contrast, the remaining solvent-rich amorphous region widens. Finally, on larger timescales, grain coarsening takes effect: surface tension contributions to the free energy are minimized, crystal interfaces suppressed, and thus percolating pathways tend to rupture.

Fig. 9 Progress of crystallization simulated for the immiscible system represented by the phase diagram Fig. 8a at T = 333 K. The volume fraction field ϕ is shown at different times for blend ratios (a) ϕ₀ = 0.4 and (b) ϕ₀ = 0.6. The simulation parameters are identical to those listed in Table 1, except M₀ = 1000 s⁻¹ and χ_aa = 1.2648. In (a), following an initial spinodal decomposition, crystals nucleate and grow during the coarsening stage of the amorphous–amorphous phase separation (top row). The crystallization is favored in the domains where the solute is in majority, here the background matrix. Subsequently, grain coarsening also takes place after crystal impingement (bottom row). In (b), thermal fluctuations trigger crystal nuclei formation in an initial metastable mixed amorphous phase. This causes AAPS by a nucleation and growth mechanism. Therefore, solvent-rich droplets form next to crystals. Nucleation and growth of both crystals and amorphous droplets proceeds until the whole mixed amorphous phase is consumed. Finally, further domain coarsening is also visible in the late stages of the process.

In comparison with Fig. 9a, Fig. 9b shows the progress of crystallization for a binary system that is metastable in the mixed amorphous state. As can be seen on the first snapshot, thermal fluctuations trigger the formation of nuclei in the initial phase which is still homogeneously mixed. Due to the associated latent heat release, the free energy level changes next to these crystal germs and the threshold for NG AAPS is locally exceeded. As a result, droplets of solvent appear at these locations. Both crystals and amorphous droplets continue to nucleate and grow until the initial homogeneous phase is fully consumed. As opposed to crystals, droplets coalesce when they impinge. For the blend ratio used in Fig. 9b, the lever rule implies a higher volume proportion of the solid phase. Thus, once the phase separation is completed, crystalline domains end up in a well-connected branched configuration which covers most of the system. Conversely, the solvent-rich amorphous regions are isolated. This still holds at later coarsening stages even though the domains become larger and certain ramifications are suppressed.

5.3 Diffusion-limited crystallization in immiscible blends

Special types of morphologies can also develop in immiscible systems where crystallization is fast in comparison to diffusion (i.e. simultaneously with relatively large M(ϕ)/Λ(ϕ)). In that case, since AAPS kinetics are directly related to the diffusion rate, crystal nucleation and growth might proceed faster than amorphous–amorphous demixing. Different formation pathways can be distinguished under these circumstances. On the one hand, at low solute volume fractions (Fig. 10a), where crystallization is only possible following a two-step nucleation route (see free energy landscape: Fig. 8b), AAPS still has to take place first. Crystals then start forming as soon as the solute content in the solute majority phase is high enough to allow for direct nucleation, even though the binodal compositions are not reached. Because the crystals grow surrounded with depletion zones due to diffusion limitations, wormlike structures that might have arisen during SD break up and the morphology evolves towards isolated, droplet-like crystals, which are ultimately subject to Ostwald ripening.

Fig. 10 Progress of crystallization for an immiscible binary blend subject to diffusion limitations at T = 333 K. The volume fraction field ϕ is shown at different times for blend ratios (a) ϕ₀ = 0.2 and (b) ϕ₀ = 0.4. For both simulations, the thermodynamic and kinetic parameters are identical to those listed in Table 1, except M₀ = 10⁶ s⁻¹ and χ_aa = 1.2648. In (a), a two-step crystallization process is undergone: the transformation is initiated by a spinodal decomposition followed by crystal nucleation and diffusion-controlled growth in the emerging solute majority phase. Conversely in (b), direct crystal nucleation and growth occurs while a parallel spinodal decomposition takes place in the remaining mixed amorphous phase. Due to diffusion limitations, the growing crystals are surrounded by depletion zones. The arising amorphous domains where solute is in majority present favorable conditions for secondary crystals to develop and the system progressively becomes more nucleation-dominated. Therefore, smaller and closer-spaced crystals form in the late stages.

On the other hand, if ϕ₀ is large enough for a one-step process to happen, crystal nucleation and growth occurs before any amorphous phase dissociation (Fig. 10b). Composition gradients within depletion zones trigger AAPS in a particular concentric pattern that is centered on crystals, as already observed in previous PF studies of placed initial crystal seeds in non-diffusion-controlled situations.¹⁶ Here, this effect is however confined to the vicinity of the crystals. The emerging high solute concentration rings that surround the depleted areas provide a preferential location for secondary germs to nucleate. Since depletion zones form as well around these new crystals, the rings with high solute content from which they originate are disrupted. Therefore, the secondary crystals still grow relatively isolated. In later stages, upon further progress of AAPS in the amorphous region, crystallization becomes more nucleation-dominated in the solute majority phase. Thus, closer and smaller-sized germs materialize with a comparatively accelerated formation rate. Finally, after the solute material is fully consumed, the global morphology exhibits a relatively wide distribution in crystalline droplet sizes. In this way, those which appeared first under more growth-dominated conditions had time to become larger, adopt rounder shapes and are well-dispersed. In contrast, crystals which emerged last are smaller, closer to each other and slightly more distorted.

5.4 Dilution-enhanced crystallization

For all simulations presented in the previous sections, the local incorporation rate of amorphous solute into crystals is assumed constant over the whole composition range (i.e. M(ϕ) = M₀). With the aim to model organic solvent–solute systems, it has nevertheless to be considered that the molecular mobility of crystalline materials such as polymers should be globally enhanced in the diluted state.²⁴ Correspondingly, a faster crystal build-up should be kinetically enabled at lower solute volume fractions. To account for this, the value of the Allen–Cahn mobility parameter M(ϕ), which locally controls the amorphous-to-crystalline conversion rate, can be indexed as follows on the self-diffusion coefficient of the solute D^(self)₁, which, in turn, depends on ϕ (see eqn (17) and SI-A, ESI†):^14,24

The transformation kinetics in simulations implementing this modification are depicted in Fig. 11. In comparison with the regular case (Fig. 4), the process does not monotonously accelerate with higher solute volume fractions, but rather exhibits a maximal crystallization rate at intermediate ϕ₀. This is explained by the presence of two driving mechanisms for crystallization, which possess opposite trends with respect to solute content. In this way, the thermodynamic driving force (eqn (3)), that becomes stronger for high solute concentrations, competes with the molecular mobility enhancement (eqn (19)) that is more effective for high dilutions (i.e. low ϕ). The exact solute volume fraction for which the crystallization rate is maximized then depends on material specific properties such as the heat of fusion or the self-diffusivities.

Fig. 11 (a) Transformation kinetics in the dilution-enhanced crystallization regime at different solvent–solute blend ratios. (b) Comparison of the crystallization half-time as a function of blend composition in the dilution-enhanced and regular crystallization regimes.

Furthermore, morphological features arising in this crystallization scenario can significantly differ from the non-dilution-enhanced case. Higher nucleation rates, which lead to smaller-sized crystals are observed for a given blend ratio. Globally, the simulated systems however still tend to be more growth-dominated towards low ϕ₀. In addition, the crystal build-up is also less homogeneous due to local variations of M with ϕ. If the system is simultaneously subject to composition-dependent M(ϕ) (eqn (19)) and diffusion-limitations, depletion zones and bead-chain like crystalline structures are obtained as well. If the mixture is immiscible, the comparatively higher nucleation rate, combined with the progressively slowing down crystal purification (as M(ϕ) decreases when the solute concentration ϕ increases), impact the interplay between crystallization and the concurrently occurring AAPS. The influence of M(ϕ) on the resulting crystal arrangements and size distributions is particularly visible when the blend is at the same time immiscible, diffusion-limited and dilution-enhanced. A summary of the corresponding morphologies is provided in Fig. 12 and in the ESI† (SI-H). Note that when the kinetic enhancement of the crystallization rate with dilution (eqn (19)) couples with these other physics, the overall process kinetics become even more sensitive to the material properties and substantial deviations from the typical transformation plots of Fig. 11 are expected. Fig. 12 also provides a final comparative overview of the morphology formation in all crystallization regimes simulated in this work.

Fig. 12 Summary of the different crystalline morphology formation regimes simulated in this study for the blend ratio ϕ₀ = 0.4 at T = 333 K. On the one hand, the rows of the table are associated with the miscibility of the mixture in the amorphous state. The system is either miscible (first row) or immiscible (second row), depending on whether the value of χ_aa is less or greater than . On the other hand, the columns refer to the limiting physical mechanism for crystal growth. Low M(ϕ)/Λ(ϕ) ratios yield crystallization scenarios where crystal growth is controlled by the crystal interface mobility (first column). In contrast, high M(ϕ)/Λ(ϕ) correspond to processes where the relatively slow diffusion of solute molecules towards crystal interfaces is limiting the crystallization (second column). Within a cell, both snapshots of the first sub-row depict the time-evolution of the crystalline structure in a blend without dilution-enhancement (M(ϕ) = M₀). In comparison, the second sub-row shows the progress of crystallization when the system is additionally dilution-enhanced (eqn (19)).

6 Discussion – comparison with experimental results for OPV systems

In order to evaluate how the present PF approach can be applied to solution-processed PALs, reports discussing the phenomenology of crystallization in material blends relevant for OPV are drawn into comparison. It is found that the simulation outcomes are in line with many experimental results. To begin with, amorphous–crystalline polymer–polymer mixtures are known to exhibit transformation kinetics similar to the sigmoids presented in Fig. 4 where the crystallization rate increases with the overall content of crystalline material.⁶⁷ It is worth mentioning that, for these particular solute–solute configurations, diffusion properties should not vary significantly with the blend ratio.

In addition, two-step nucleation processes where crystals appear within previously formed amorphous solute-rich clusters have been theorized and observed in protein, as well as in other organic, inorganic, and colloidal solutions,^68,69 which agrees with predictions from the current PF model. The accelerated crystallization kinetics caused by an initial amorphous demixing are also consistent with publications of experiments on polymer–polymer blends.^70–72 It has been suggested that the rate of transformation is enhanced due to so-called up-hill diffusion of polymer chains into their corresponding majority phase where the thermodynamic driving force for crystallization is stronger. This corresponds to the mechanism observed in the present simulations.

Note again that for these systems the diffusion coefficients are expected to be comparable throughout the whole volume fraction range. Evidence for dilution-enhanced crystal growth rates has also been reported for polymer solutions.⁷³ However, it is specified that the dependency of the crystallization kinetics on composition can vary substantially with the investigated material system, as the change with ϕ of other blend properties (e.g. miscibility, interdiffusivity, glass transition temperature, crystal growth dimensionality, etc.) is strongly determined by the selected materials. Thereby, no general composition–crystallization rate relationship could be derived. This is in agreement with the fact that the progress of crystallization can be significantly altered when molecular mobility enhancements due to dilution are combined with diffusion-control or AAPS physics, for instance, as it has been witnessed in this work.

Regarding particular applications for OPV, mixtures involving isotropically crystallizing species, such as PCBM in the PCE11-PCBM system that was visualized by Levitsky et al. using optical micrographs and vapor phase infiltration enhanced scanning electron microscopy,⁷⁴ or the MDMO-PPV:PCBM blends imaged by Hoppe and co-workers with atomic force microscopy and scanning electron microscopy,^75–78 feature morphological structures that are similar to those reproduced in this work and provide encouraging perspectives for further investigations with the PF model as it is presented in this study. As mentioned before, anisotropic crystal growth, which is expected for a vast class of organic semi-conducting materials,⁷⁹ can be accounted for as well with this type of modelling framework.³²

Furthermore, the existence of depleted areas around crystals was evidenced with energy-filtered transmission electron microscopy for solvent–vapor annealed DRCN5T:PCBM all small molecule PALs by Harreiß et al.⁸⁰ Analyses with the current PF method can typically help in this case to determine whether diffusion-limited crystallization is responsible for the formation of such structures. Diffusion properties in the P3HT-PCBM model system leading to depletion zones around PCBM crystals were also studied by Berriman et al.⁸¹ Moreover, Yu et al.⁸² commented on the relevance of crystallization in strongly diffusion-controlled configurations for OPV systems. The importance of the changes in crystalline morphology encountered as the crystal growth rate becomes faster than the molecular diffusion rate has additionally been discussed by Wang et al.²⁰ alongside with corresponding post-treatment strategies to adopt in order to optimize the resulting PALs for PCE. Insights from simulations can be of use here to identify when diffusion limitations become detrimental (or beneficial) for device performance, so that appropriate processing measures are undertaken.

It is recognized as well that the interplay of AAPS and crystallization is crucial for BHJ manufacturing.²⁰ In general, OPV systems are likely to enter the miscibility gap of the phase diagram as composition and temperature significantly change during processing.^8,9 Thus, several studies pointed out the significance of tailoring the components' miscibility for amorphous–crystalline morphology optimization,⁵⁷ the effects of spinodal decomposition on subsequent crystallization during thermal annealing,⁷⁴ or issues with non-fullerene acceptor phase separation and pre-aggregation in solution.^37,83 As demonstrated in Section 5.2, complementary numerical investigations can be carried out with the present approach to gain a more detailed understanding of the competition between the physical phenomena involved in these processes.

Finally, it can be remarked that the computations shown here are two-dimensional, so as to reduce computational effort and allow for larger sets of simulations to be performed. The highlighted sensitivities to thermodynamic and kinetic parameters are nonetheless expected to be qualitatively similar in three dimensions. Likewise, the witnessed crystallization regimes and their respective morphological features (e.g. depletion zones, spinodal decomposition pattern, etc.) are also anticipated to arise in that case. Nevertheless, the presence of an additional spatial dimension may lead to new topological configurations. Regarding for instance bicontinuity of demixed phases, which is an important descriptor for performance of organic photoactive layers,⁸⁴ it can be recognized that more structures are likely to develop with percolating pathways because the third dimension permits more possibilities for domains to circumvent each other without breaking connectivity. Further simulations in 3D are therefore of interest to determine additional material parameter combinations which can yield useful morphologies for OPV devices.

7 Conclusions

In conclusion, the crystallization behavior in binary, non-evaporating, amorphous–crystalline mixtures was assessed with the help of numerical simulations. More specifically, the objective was to gain a detailed comprehension of possible crystalline structure formation mechanisms based on the governing physics captured by the presented Phase-Field framework for solution-processed photoactive layer simulation. In this way, several mixture properties such as the blend ratio, the strength of thermodynamic driving force for crystallization, the height of the nucleation barrier, the surface tension, the miscibility of the mixed species, or their inter- and self-diffusivities, were found to play a crucial role in the investigated morphology formation process. Upon exploration of the thermodynamic and kinetic parameter space, the sensitivity of the transformation kinetics was evaluated and various crystallization scenarios were achieved. A more exhaustive analysis of the realized diffusion-limited, demixing-assisted and dilution-enhanced crystallization regimes revealed distinct morphology formation pathways with specific process kinetics, crystal size distributions, spatial arrangements and nucleation/growth balances. Thereby, remarkable morphological features such as depletion zones forming around crystals or percolating, bicontinous domains were produced as well. Moreover, the interplay of physical phenomena, which would separately lead to different crystallization modes, was also observed to result in unique geometrical configurations. Finally, experimental results reported in the literature were drawn into comparison and found to be in good qualitative agreement with numerical simulation outcomes.

Thus, the employed phase-field modelling approach demonstrates highly promising perspectives for understanding the interaction between complex nanoscale processes that take place simultaneously during photoactive film fabrication. It is expected that theoretical insights provided by simulations will substantially help the identification of process–structure relationships and support the derivation of guidelines for optimal manufacturing conditions. In order to increase even further the prediction range for morphological features that might arise during the solution-deposition process, supplementary physical phenomena are anticipated to be of importance. Future code developments therefore include extensions for anisotropy, semi-crystallinity, substrate interactions and chemical reactions. In addition, upcoming investigations will concentrate on the validation of model predictions against experimental measurements for specific state-of-the-art material systems. To this end, particular attention will be dedicated to accurate input parameter acquisition, in situ morphological property inspection with adequate thin-film characterization techniques, and advanced quantitative comparisons using appropriate morphology descriptors. This will ultimately allow for reliable simulations of the complete photoactive film formation upon drying, which involves concurrent crystallization, demixing, evaporation and hydrodynamic processes.

Author contributions

M. Siber: methodology, investigation, formal analysis, data curation, writing – original draft, writing – review and editing. O. J. J. Ronsin: conceptualization, software, supervision, writing – review and editing. J. Harting: project administration, funding acquisition, writing – review and editing.

Data availability

In compliance with the regulations for projects funded by the German Research Foundation (DFG), the simulation data used for this article is made publicly accessible (see DOI https://doi.org/10.5281/zenodo.8422262).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge financial support by the German Research Foundation (DFG, Project HA 4382/14-1), the European Commission (H2020 Program, Project 101008701/EMERGE), and the Helmholtz Association (SolarTAP Innovation Platform).

Notes and references

1. Z. Li, T. R. Klein, D. H. Kim, M. Yang, J. J. Berry, M. F. A. M. van Hest and K. Zhu, Nat. Rev. Mater., 2018, 3, 1–20.
2. C. J. Brabec, A. Distler, X. Du, H.-J. Egelhaaf, J. Hauch, T. Heumueller and N. Li, Adv. Energy Mater., 2020, 10, 2001864.
3. F. Yang, D. Jang, L. Dong, S. Qiu, A. Distler, N. Li, C. J. Brabec and H.-J. Egelhaaf, Adv. Energy Mater., 2021, 11, 2101973.
4. R. R. Søndergaard, M. Hösel and F. C. Krebs, J. Polym. Sci., Part B: Polym. Phys., 2013, 51, 16–34.
5. S. T. Williams, A. Rajagopal, C.-C. Chueh and A. K.-Y. Jen, J. Phys. Chem. Lett., 2016, 7, 811–819.
6. Y. Wang, J. Lee, X. Hou, C. Labanti, J. Yan, E. Mazzolini, A. Parhar, J. Nelson, J.-S. Kim and Z. Li, Adv. Energy Mater., 2021, 11, 2003002.
7. N. E. Jackson, J. Phys. Chem. B, 2021, 125, 485–496.
8. C. McDowell, M. Abdelsamie, M. F. Toney and G. C. Bazan, Adv. Mater., 2018, 30, 1707114.
9. Z. Peng, N. Stingelin, H. Ade and J. J. Michels, Nat. Rev. Mater., 2023, 1–17.
10. O. Wodo and B. Ganapathysubramanian, Comput. Mater. Sci., 2012, 55, 113–126.
11. S. Kouijzer, J. J. Michels, M. van den Berg, V. S. Gevaerts, M. Turbiez, M. M. Wienk and R. A. J. Janssen, J. Am. Chem. Soc., 2013, 135, 12057–12067.
12. C. Schaefer, J. J. Michels and P. van der Schoot, Macromolecules, 2016, 49, 6858–6870.
13. J. J. Michels, K. Zhang, P. Wucher, P. M. Beaujuge, W. Pisula and T. Marszalek, Nat. Mater., 2021, 20, 68–75.
14. O. J. J. Ronsin and J. Harting, Adv. Theory Simul., 2022, 2200286.
15. B. König, O. J. J. Ronsin and J. Harting, Phys. Chem. Chem. Phys., 2021, 23, 24823–24833.
16. O. J. J. Ronsin and J. Harting, Energy Technol., 2020, 8, 1901468.
17. O. J. J. Ronsin, D. Jang, H.-J. Egelhaaf, C. J. Brabec and J. Harting, Phys. Chem. Chem. Phys., 2020, 22, 6638–6652.
18. O. J. J. Ronsin, D. Jang, H.-J. Egelhaaf, C. J. Brabec and J. Harting, ACS Appl. Mater. Interfaces, 2021, 13, 55988–56003.
19. O. J. J. Ronsin and J. Harting, ACS Appl. Mater. Interfaces, 2022, 14, 49785–49800.
20. Z. Wang, K. Gao, Y. Kan, M. Zhang, C. Qiu, L. Zhu, Z. Zhao, X. Peng, W. Feng, Z. Qian, X. Gu, A. K.-Y. Jen, B. Z. Tang, Y. Cao, Y. Zhang and F. Liu, Nat. Commun., 2021, 12, 332.
21. L. Gránásy, G. I. Tóth, J. A. Warren, F. Podmaniczky, G. Tegze, L. Rátkai and T. Pusztai, Prog. Mater. Sci., 2019, 106, 100569.
22. O. Wodo and B. Ganapathysubramanian, J. Comput. Phys., 2011, 230, 6037–6060.
23. D. Zhou, A.-C. Shi and P. Zhang, J. Chem. Phys., 2008, 129, 154901.
24. C.-S. Kim, D. M. Saylor, M. K. McDermott, D. V. Patwardhan and J. A. Warren, J. Biomed. Mater. Res., Part B: Appl. Biomater., 2009, 90B, 688–699.
25. D. M. Saylor, J. E. Guyer, D. Wheeler and J. A. Warren, Acta Biomater., 2011, 7, 604–613.
26. D. M. Saylor, C. Forrey, C.-S. Kim and J. A. Warren, Ann. Biomed. Eng., 2016, 44, 548–559.
27. R. A. Matkar and T. Kyu, J. Phys. Chem. B, 2006, 110, 12728–12732.
28. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
29. P. C. Hohenberg and A. P. Krekhov, Phys. Rep., 2015, 572, 1–42.
30. J. I. Lauritzen and J. D. Hoffman, J. Res. Natl. Bur. Stand., Sect. A, 1960, 64A, 73–102.
31. T. Takaki, ISIJ Int., 2014, 54, 437–444.
32. L. Gránásy, L. Rátkai, A. Szállás, B. Korbuly, G. I. Tóth, L. Környei and T. Pusztai, Metall. Mater. Trans. A, 2014, 45, 1694–1719.
33. P. J. Flory, Principles of Polymer Chemistry, Cornell University Press, 1953.
34. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, Oxford, New York, 2003.
35. J. Perea, S. Langner, M. Salvador, B. Sanchez, N. Li, C. Zhang, G. Járvás, J. Kontos, A. Dallos, A. Aspuru-Guzik and C. Brabec, J. Phys. Chem. C, 2017, 121, 18153–18161.
36. L. Ye, H. Hu, M. Ghasemi, T. Wang, B. A. Collins, J.-H. Kim, K. Jiang, J. H. Carpenter, H. Li, Z. Li, T. McAfee, J. Zhao, X. Chen, J. L. Y. Lai, T. Ma, J.-L. Bredas, H. Yan and H. Ade, Nat. Mater., 2018, 17, 253–260.
37. G. Yao, Y. Ge, X. Xiao, L. Zhang, N. Yi, H. Luo, S. Yuan and W. Zhou, ACS Appl. Energy Mater., 2022, 5, 1193–1204.
38. W. Hu, Polymer Physics, Springer Vienna, Vienna, 2013.
39. S. M. Allen and J. W. Cahn, Acta Metall., 1979, 27, 1085–1095.
40. C. Shen, J. Simmons and Y. Wang, Acta Mater., 2007, 55, 1457–1466.
41. L. Gránásy, T. Pusztai, G. Tegze, J. A. Warren and J. F. Douglas, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011605.
42. D. Lee, J.-Y. Huh, D. Jeong, J. Shin, A. Yun and J. Kim, Comput. Mater. Sci., 2014, 81, 216–225.
43. S. R. D. Groot and P. Mazur, Non-Equilibrium Thermodynamics, Courier Corporation, 1984.
44. M. Hillert, Phase Equilibria, Phase Diagrams and Phase Transformations: Their Thermodynamic Basis, Cambridge University Press, Cambridge, UK; New York, 2008.
45. J. W. Cahn, J. Chem. Phys., 1965, 42, 93–99.
46. L. Onsager, Phys. Rev., 1931, 37, 405–426.
47. E. J. Kramer, P. Green and C. J. Palmstrøm, Polymer, 1984, 25, 473–480.
48. P. G. de Gennes, J. Chem. Phys., 1980, 72, 4756–4763.
49. A. Z. Akcasu, Macromol. Theory Simul., 1997, 6, 679–702.
50. L. Wolff, S. H. Jamali, T. M. Becker, O. A. Moultos, T. J. H. Vlugt and A. Bardow, Ind. Eng. Chem. Res., 2018, 57, 14784–14794.
51. H. Cook, Acta Metall., 1970, 18, 297–306.
52. R. Horst, Macromol. Theory Simul., 1995, 4, 449–458.
53. R. Horst, Macromol. Theory Simul., 1996, 5, 789–800.
54. D. R. Kozub, K. Vakhshouri, L. M. Orme, C. Wang, A. Hexemer and E. D. Gomez, Macromolecules, 2011, 44, 5722–5726.
55. Y. Zhu, A. Gadisa, Z. Peng, M. Ghasemi, L. Ye, Z. Xu, S. Zhao and H. Ade, Adv. Energy Mater., 2019, 9, 1900376.
56. R. A. Matkar and T. Kyu, J. Phys. Chem. B, 2006, 110, 16059–16065.
57. L. Ye, B. A. Collins, X. Jiao, J. Zhao, H. Yan and H. Ade, Adv. Energy Mater., 2018, 8, 1703058.
58. M. Ghasemi, H. Hu, Z. Peng, J. J. Rech, I. Angunawela, J. H. Carpenter, S. J. Stuard, A. Wadsworth, I. McCulloch, W. You and H. Ade, Joule, 2019, 3, 1328–1348.
59. A. Wadsworth, Z. Hamid, J. Kosco, N. Gasparini and I. McCulloch, Adv. Mater., 2020, 32, 2001763.
60. B. Schmidt-Hansberg, M. Sanyal, M. F. G. Klein, M. Pfaff, N. Schnabel, S. Jaiser, A. Vorobiev, E. Müller, A. Colsmann, P. Scharfer, D. Gerthsen, U. Lemmer, E. Barrena and W. Schabel, ACS Nano, 2011, 5, 8579–8590.
61. S. Ludwigs, P3HT Revisited – From Molecular Scale to Solar Cell Devices, Springer Berlin Heidelberg, Berlin, Heidelberg, 2014, vol. 265.
62. P. R. Berger and M. Kim, J. Renewable Sustainable Energy, 2018, 10, 013508.
63. Q. Guo, Q. Guo, Y. Geng, A. Tang, M. Zhang, M. Du, X. Sun and E. Zhou, Mater. Chem. Front., 2021, 5, 3257–3280.
64. M. Avrami, J. Chem. Phys., 1939, 7, 1103–1112.
65. M. Avrami, J. Chem. Phys., 1940, 8, 212–224.
66. C. A. Kelly and M. J. Jenkins, Polym. J., 2022, 54, 249–257.
67. L. Mandelkern, Crystallization of Polymers: Kinetics and Mechanisms, Cambridge University Press, Cambridge, 2nd edn, 2004, vol. 2, pp. 282–347.
68. D. Erdemir, A. Y. Lee and A. S. Myerson, Acc. Chem. Res., 2009, 42, 621–629.
69. P. G. Vekilov, Nanoscale, 2010, 2, 2346–2357.
70. M. Tsuburaya and H. Saito, Polymer, 2004, 45, 1027–1032.
71. K. Sugeno, S. Kokubun and H. Saito, Polymers, 2020, 12, 2730.
72. K. Takamatsu and H. Saito, Polymer, 2022, 238, 124418.
73. L. Mandelkern, Crystallization of Polymers: Kinetics and Mechanisms, Cambridge University Press, Cambridge, 2nd edn, 2004, vol. 2, pp. 380–432.
74. A. Levitsky, S. A. Schneider, E. Rabkin, M. F. Toney and G. L. Frey, Mater. Horiz., 2021, 8, 1272–1285.
75. H. Hoppe, M. Niggemann, C. Winder, J. Kraut, R. Hiesgen, A. Hinsch, D. Meissner and N. S. Sariciftci, Adv. Funct. Mater., 2004, 14, 1005–1011.
76. H. Hoppe, M. Drees, W. Schwinger, F. Schaffler and N. Sariciftcia, Synth. Met., 2005, 152, 117–120.
77. H. Hoppe and N. S. Sariciftci, J. Mater. Chem., 2006, 16, 45–61.
78. P. A. Troshin, H. Hoppe, A. S. Peregudov, M. Egginger, S. Shokhovets, G. Gobsch, N. S. Sariciftci and V. F. Razumov, ChemSusChem, 2011, 4, 119–124.
79. C. Wang, H. Dong, L. Jiang and W. Hu, Chem. Soc. Rev., 2018, 47, 422–500.
80. C. Harreiß, S. Langner, M. Wu, M. Berlinghof, S. Rechberger, J. Will, M. Conroy, U. Bangert, T. Unruh, C. J. Brabec and E. Spiecker, Solar RRL, 2022, 2200127.
81. G. A. Berriman, J. L. Holdsworth, X. Zhou, W. J. Belcher and P. C. Dastoor, AIP Adv., 2015, 5, 097220.
82. L. Yu, D. Qian, S. Marina, F. A. A. Nugroho, A. Sharma, S. Hultmark, A. I. Hofmann, R. Kroon, J. Benduhn, D.-M. Smilgies, K. Vandewal, M. R. Andersson, C. Langhammer, J. Martín, F. Gao and C. Müller, ACS Appl. Mater. Interfaces, 2019, 11, 21766–21774.
83. J. Xue, H. Bilal Naveed, H. Zhao, B. Lin, Y. Wang, Q. Zhu, B. Wu, Z. Bi, X. Zhou, C. Zhao, K. Zhou and W. Ma, J. Mater. Chem. A, 2022, 10, 13439–13447.
84. O. Wodo, J. Zola, B. S. Sarath Pokuri, P. Du and B. Ganapathysubramanian, Mater. Discovery, 2015, 1, 21–28.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3tc03047d

---