Phase Separation in PCDTBT:PCBM Blends: from Flory-Huggins Interaction Parameters to Ternary Phase Diagrams

The substantial increase in the efficiency of organic solar cells achieved in recent years would not have been possible without work on the synthesis of new materials and understanding the relationship between the morphology and performance of organic photovoltaic devices. The structure of solvent-cast active layers is a result of phase separation in mixtures of donor and acceptor components. To a large extent, this process depends on the interactions between the components of the mixture. Here, we present a systematic analysis of the morphology of poly[N-9′-heptadecanyl-2,7-carbazole-alt-5,5-(4′,7′-di-2-thienyl-2′,1′,3′-benzothiadiazole)] (PCDTBT) and [6,6]-phenyl-C71-butyricacid methyl ester (PC70BM) films in terms of the ternary phase diagram. The interaction parameters between PCDTBT and four different solvents, namely chloroform, chlorobenzene, o-dichlorobenzene, and toluene, were estimated based on swelling experiments. Based on these values, ternary phase diagrams of PCDTBT:PC70BM in different solvents were calculated. The morphology of spin-coated films with different blend ratios cast from different solvents is discussed in terms of the obtained phase diagrams.


INTRODUCTION
Polymer blends have already found applications in numerous fields, such as organic electronics, optics, and biotechnology, [1] mainly due to their ability to form a variety of structures, e.g., lamellar, lateral, regular, or hierarchic, [2] with a large size range from nano-to micrometers. Additionally, compared to metal or classic semiconductors, polymer blends can be dissolved and deposited as thin films in a one-step procedure by a coating, printing, or roll-to-roll technique. During solvent evaporation, the polymer blends undergo phase separation, resulting in the formation of different domains. Their final size and shape depend on many thermodynamic (temperature, pressure), processing (coating speed) and material parameters (solution concentration, solubility parameter, miscibility of components, parameter). [1,3] Because the performance of devices based on polymer blends is sensitive to the final phase domain structure, the key issue in device optimization is to understand and control the phase separation process. In the case of organic photovoltaics (OPV), the optimal phase domain size of the active layer is approximately 10 nm and should resemble a comb-like structure. This restriction is linked to a mechanism in which free charges are formed. In contrast to their inorganic counterparts, light absorption leads first to creation of Frenkel exciton (electron-hole pair bounded with Coulombic attraction 0.3-0.5 V [4] ). To split the excitons into the free electrons and holes, it is necessary to overcome the binding energy, which occurs when the exciton reaches the donor-acceptor interface. However, the exciton diffusion path amounts to only 10 nm. [5] Therefore, to maximize the charge extraction from the OPV, a donor-acceptor interpenetrating network should be prepared. In such a system, the exciton can reach the donoracceptor interface within its lifetime. Additionally, the interpenetrating network should form a continuous pathway to ensure that the free charges can reach their respective electrodes. In many organic solar cells, semiconducting polymers absorbing in the visible light spectrum are used as donors, while fullerene derivatives, such as [6,6]-phenyl-C 61 -butyric acid methyl ester (PC60BM) or [6,6]-phenyl-C 71 -butyric acid methyl ester (PC70BM), are used as acceptors. Although the last decades witnessed many improvements to OPV morphology through experimental methods, [6−8] several recent studies have been dedicated to modeling and predicting the polymer:fullerene blend behavior, [9−15] typically focusing on specific donor/acceptor/solvent ternary blends. Such research is usually based on an analysis of the Gibbs ternary phase diagrams.
The starting point to construct a ternary phase diagram is to determine the Flory-Huggins interaction parameter χ between the polymer and solvents. Typically, the interaction parameter χ, or more precisely its purely enthalpic concentration independent fraction is estimated based on the solubility parameters δ determined experimentally via osmosis, vapor pressure measurement, inverse gas chromatography, or contact angle. [16] Light scattering experiments allow to determine the interaction parameter between polymer and solvent directly in the limit of low polymer concentration. [17] Recently, it was shown that the interaction parameters of semiconducting polymers in the high polymer concentration limit can be estimated from swelling experiments. [18] A polymer layer swells in the presence of solvent vapors because solvent molecules diffuse into the polymer film. The extent of swelling, according to the regular solution, [19] depends directly on the value of the interaction parameter between polymers and solvents.
A substantial increase in the efficiency of organic solar cells has been achieved in recent years and devices with power conversion efficiency (PCE) higher than 14% were reported. [20,21] This would not have been possible without the enormous progress in the synthesis of new acceptors and low bandgap polymeric semiconductors. [22] One of the new commercially available low bandgap polymers, poly[N-9′-heptadecanyl-2,7-carbazole-alt-5,5-(4′,7′-di-2-thienyl-2′,1′,3′-benzothiadiazole)], PCDTBT, has aroused great interest as a donor material for bulk heterojunction solar cells. The solar cells based on the PCDTBT:PC70BM mixture are characterized by a high open-circuit voltage, high efficiency up to 7.5%, and a possible life time up to 7 years. [23] Moreover, Park et al. [24] showed that the internal quantum efficiency can be close to 100%, which means that almost all absorbed photons generate carriers, which are subsequently collected by the electrodes. Despite the fact that PCDTBT is widely used, little attention has been paid to quantifying an interaction parameter between PCDTBT and solvents. Information on this parameter can help to predict the phase separation behavior in donor-acceptor mixtures and streamline the OPV optimization process.
Here, we present an analysis of a (PCDTBT):fullerene (PC70BM) blend behavior via ternary phase diagrams. The interaction parameters between PCDTBT and four different solvents, namely chloroform, chlorobenzene, o-dichlorobenzene, and toluene, were estimated based on swelling experiments. The obtained values were used to calculate spinodal and binodal lines limiting the unstable and one-phase region in the blends of PCDTBT:PC70BM dissolved in the four chosen solvents. Finally, PCDTBT:PC70BM films with different blend ratios in different solvents were prepared by spin-coating. Their morphology was discussed based on analysis of the ternary phase diagrams. In this experiment, we concentrated on a particular blend system, PCDTBT and PC70BM dissolved in several solvents, which is currently one of the most widely studied donor-acceptor mixtures used in the field of bulk heterojunction solar cells. [25] Despite the popularity of PCDTBT, to our knowledge, this is the first time PCDTBT interaction parameters with different solvents were measured and respective PCDTBT:PC70BM:solvent mixtures were analyzed via a ternary phase diagram.

Sample Preparation
Neat PCDTBT polymer films were spin-cast from chlorobenzene on top of silicon substrates. Silicon dies ca. 1 cm × 1 cm with native silicon oxide were used as substrates for samples dedicated to determining the refractive index n and extinction coefficient k with spectral ellipsometry. Dies with thermally grown SiO 2 were used as substrates in swelling experiments. All substrates were cleaned by sonication in toluene.

Swelling of Polymer Films
A laboratory-built system composed of two gas flow controllers (MKS), and a white light reflectance spectrometer (WLRF) was used to study the swelling of polymer films exposed to vapor of different volatile solvents. The fraction of solvent vapor in the atmosphere above the sample was adjusted by mixing nitrogen flux passing a Drechsler's bottle filled with volatile solvent and pure N 2 . Variations in film thickness were recorded in real time using a FRBasic spectrometer (Thetametrisis, Athens Greece). The refractive index n and extinction coefficient k for PCDTBT films as a function of light wavelength were determined prior to the swelling experiment using spectroscopic ellipsometry (Senetech model SE800).

Morphology
An Agilent 5500 atomic force microscope was used to study the morphology of the PCDTBT:PC70BM blend films. Topographic images with dimensions of 25 μm × 25 μm and 50 μm × 50 μm were acquired in contact/intermittent contact mode. Proportional and integral gains were adjusted individually for each scan to obtain pictures with the best quality.

Swelling of PCDTBT
A schematic of the laboratory-built apparatus for the swelling experiment is presented in Fig. 1. The polymer film under examination was placed in a chamber flushed with nitrogen carrying a predetermined fraction of solvent vapor. The composition of the atmosphere over the sample was adjusted by mixing gas flux bubbled through the chosen solvent filling Drechsler's bottle with a flux of pure nitrogen. The ratio of the fluxes was adjusted by two gas flow controllers while the total gas flow was kept constant at 2000 mL/min; the gas flow was adjusted experimentally to have neglectable effect on the final results of the swelling experiments. The expansion of the polymer films under solvent uptake was traced by a white light reflectance interferometer. Fringes in the spectrum of the white light reflected from the PCDTBT cast on top of silicon with thermally grown SiO 2 were recorded at 5-s intervals. The interference spectra were compared with the reference spectrum of the halogen lamp to yield, in real time, the optical path in the polymer film. The refractive index, n, and extinction coefficient, k, determined for PCDTBT (see the electronic supplementary information, ESI) were used to calculate the relative expansion of film thickness d/d 0 . The initial thickness d 0 was calculated as an average of the thicknesses measured in the first 10 min of experiment when the solvent vapor partial pressure was close to zero.
Polymer films were exposed to the sequence of N 2 fluxes carrying the solvent with different relative partial pressures p/p sat inside the measurement chamber. During the experiment, the relative partial pressure in the chamber varied from 0.05 to 0.85. Each sequence was separated by pure N 2 fluxes when the polymer film relaxed to its original thickness. Penetration of solvent molecules into the polymer film resulted in an increase in the distance between polymer chains [26] and expansion of the film. For a long exposure of the polymer layer to the solvent vapor, the system should equilibrate, and the film thickness should stabilize at value d eq . Fig. 2 presents the whole sequence of relative expansions d/d 0 of the PCDTBT layer exposed to volatile components with different relative partial pressures in the chamber. Qualitatively, the film expansion is related to the compatibility between polymer and solvent; for more compatible solvent, higher expansion is observed for the same relative partial pressure. For each solvent under examination, a rapid increase in the relative thickness of the polymer film was ob-

Fig. 2
Relative expansion of PCDTBT films exposed to the sequence of various solvent vapors: chlorobenzene, chloroform, o-dichlorobenzene, and toluene with different relative partial pressures.
served, but the swelling dynamic depended on the type of diluent. In the case of chloroform and toluene, for higher relative pressure, after swelling rapidly, the layer started to deflate. After longer observation time, the regression of the swelling in the toluene and chloroform vapors could be explained by a drop in the saturated vapor pressure carried by the nitrogen flux. Rapid evaporation of more volatile liquids (chloroform, toluene) under nitrogen bubbling led to a visible lowering of the fluid level in the Drechsler's bottle and its cooling. For less volatile solvents (chlorobenzene, dichlorobenzene), those effects were unnoticeable, and an ongoing increase in thickness was observed for each exposure condition. To calculate the film thickness at equilibrium with the vapor for relative partial pressures in the chamber, the following model was used: [27] where d eq is the film thickness in equilibrium, and A and α are characteristic constants. Fig. S2(a) (in ESI) shows that the experimental data obtained for low relative pressure (p/p sat < 0.45) reflect the exponential convergence (1). For high relative pressure, only the initial part of the data, corresponding to an increasing thickness (Fig. 2), was fitted with Eq. (1).

Flory-Huggins Interaction Parameter and Solubility Parameter
Sorption of solvent molecules by nonpolar or weakly polar polymers can be regarded as a dissolution process analogous to that in liquid solvents and can be described analogously to the regular solution theory: [19,28] ln ( relating the difference in chemical potential of VC vapor over swollen polymer (polymer solution) μ VC and over a pure VC liquid with relative partial pressure p/p sat . In equilibrium between the film solid and volatile component, the polymer volume fraction is given by the reverse of relative film expansion Φ P = d 0 /d eq , and taking into account that molar volumes of volatile compounds V VC are much smaller than the polymer molar volume V P , Eq. (2) can be approximated as: [18] ln ( The values of relative film expansion d eq /d 0 estimated from swelling experiments were used to calculate the Flory-Huggins interaction parameters between PCDTBT and the selected solvents according to Eq. (3). The results are presented in Fig. 3.
In the mean field model, interactions between mixture components, such as polymer and solvent, depend on the solubility parameters δ i that characterize each of them (i = VC, P for solvent and polymer, respectively). The solubility parameter is simply the square root of the cohesive energy density. This approach was originally developed to predict the mixing of simple nonpolar solvents and was later extended to polar solvents and polymers. [16] The most direct way to estimate δ is to measure the energy of vaporization, but this is possible only for volatile materials like solvents. In the case of polymers, indirect methods must be used, such as solvent testing, osmotic pressure measurement, surface tension measurement, or inverse phase gas chromatography (for more details, see Ref. [16]). Additionally, there are some different computational tools [29] based on a group contribution technique, which assumes that the total energy of vaporization of a polymer is a sum of the energy of vaporization of different functional groups. Here, we estimated the solubility parameter of PCDTBT by taking into account the interaction parameters obtained in swelling experiments and the relation: [30] The solubility parameter and molar volume of the solvent were taken from Ref.

Ternary Phase Diagrams and Phase Separation in PCDTB:PC70BM
The Flory-Huggins mean field theory [33] assumes that the enthalpic contribution to the free energy of mixing of a system composed of two types of macromolecules (A and B) dissolved in common solvent (VC) depends on the mutual interactions between the components described by interaction parameters χ ij (where i ≠ j = A, B, or VC). The Gibbs free energy of mixing is given by: 5) where N A and N B are the number of sites occupied by macromolecules of type A and B, respectively. Following the idea given by Nilsson et al., [12] we calculated ternary phase diagrams to predict the phase behavior of the PCDTBT:PC70BM blend in four different solvents: chloroform, chlorobenzene, odichlorobenzene, and toluene. The interaction parameters were calculated based on the solubility parameter of PCDTBT estimated here and known from the literature for the other blend components; details are given in ESI. All the numeric values used in the calculations are summarized in Table 1.
The ternary phase diagrams of the PCDTBT:PC70BM blend in four different solvents are presented in Fig. 5. The blue line limiting the unstable regime on each Gibbs triangle, the socalled spinodal line, connects points for which the second derivative of the Gibbs free energy is equal to zero. Inside the unstable region, the mixture decomposes spontaneously. The coexistence curves (binodals) limiting the stable regions of the phase diagrams are plotted in red. To calculate binodal compositions, we developed an algorithm described in detail in Appendix I. In brief, the algorithm is based on a numerical solution of a set of nonlinear differential equations (see A4). The calculation starts at the critical point where the binodal is tangent to the spinodal line, and iteratively, two branches of the coexistence curve are calculated.
The phase diagrams of the PCDTBT:PC70BM mixture dissolved in chloroform and chlorobenzene are similar (Figs. 5a and 5b), with the critical point shifted off the center of the Gibbs triangle. Both binodal and spinodal lines are asymmet- Table 1 Flory-Huggins interaction parameter χ ij for ternary blend solvent (1), PCDTBT (2), and PC70BM (3). In each case, the effective degree of polymerization amounts was N 1 = 1, N 2 = 440, and N 3 = 6.  Table 1. The red lines represent the binodal and blue lines the spinodal composition. Additionally, starting at the right corner, the direction of the solvent quench for the 2:1, 1:1, 1:3, and 1:4 PCDTBT:PC70BM ratios is marked.
ric, which is a result of the asymmetry in the molecular size of the polymer and fullerene derivative. The one-phase region is restricted to mixtures enriched in PCDTBT and with a low content of PC70BM. The critical point of the blend dissolved in odichlorobenzene is shifted towards more dilute solutions. The single-phase region visible in Fig. 5(c) is comparatively smaller than for chloroform or chlorobenzene solutions, and for the blend dissolved in toluene, no one-phase region is observed (Fig. 5d). For chloroform and chlorobenzene, the interaction parameter between solvent and polymer χ 12 is relatively low (0.431 for CF and 0.178 for CB), indicating that both can be treated as good solvents for PCDTBT. The shift in the critical point in o-dichlorobenzene solution is related to the increase in the interactions between the polymer and solvent for χ 12 = 0.544. The lack of a single-phase area in Fig. 5(d) results from strong repulsive interactions between PCDTBT and toluene and PC70BM and toluene, where χ 12 = 0.942 and χ 13 = 0.929, respectively. Optospinograms presented in Fig. S4 (in ESI) suggest that the area of stability of the real PCDTBT: PCB70BM:toluene mixture is indeed limited to mixtures with high solvent concentration. In many technological applications, thin polymer films are formed by solvent casting. First, a wet film is formed by spinor blade-coating, and it solidifies as a result of solvent evaporation. In the case of a ternary blend (macromolecules of two types and solvents), the average composition of the drying film can be represented by a straight line connecting the corner of the Gibbs triangle representing the solvent with a point on opposite side corresponding to the composition of the dissolved compounds. When the solvent evaporates, the whole system moves towards the higher concentration along the respective straight line until it reaches the border of the unstable region. After the border, the mixture is no longer homogeneous, and demixing occurs via spinodal decomposition or nucleation-and-growth, depending on the location of the intersection point. As noted by Nilsson and coworkers, [12] the position at which the system crosses the spinodal line on the phase diagram influences the final morphology of vitrifying films.
The dashed lines starting from the right corner of each triangle presented in Fig. 5 indicate the direction of the PCDTBT:PC70BM solvent quenched with different blend weight ratios, namely 2:1, 1:1, 1:3, and 1:4. The topography of films formed by spin-casting under conditions represented by these lines is presented in Fig. 6, and corresponding depth profiles showing composition variation in the direction perpendicular to the free surface are presented in Fig. S5 (in ESI). Both the polymer:fullerene mixing ratios and the type of solvent affect the final morphology. Phase separation in the PCDTBT:PC70BM blend cast from different solvents leads to a wide range of features: distinct islands for toluene, lateral structures for chloroform, nanoscale separation for chloro- benzene, and a homogeneous layer for o-dichlorobenzene, chlorobenzene and toluene.
The largest variety of structures is observed for the toluene cast mixture. For the mixing ratio PCDTBT:PC70BM 2:1, a nearly flat surface with very fine features is observed (roughness, measured as root mean square (RMS) ≈ 0.44 nm), and respective depth profiles reveal homogeneous distribution of both components in continuous films. With increasing PC70BM content, the islands with circular shapes become more pronounced; as seen in corresponding profiles, the silicon signal raises already at the beginning of sputtering process indicating that the polymer layers are not continuous. For a ratio of 1:3, the islands are 66 ± 2 nm high and 543 ± 10 nm wide, and for 1:4, they are even higher and wider at 127 ± 3 nm and 848 ± 34 nm, respectively. Similar islands of PCBM were observed by Hoppe et al. [34] for a MDMO-PPV:PC60BM film spin-coated from a toluene solution. Their size was one order of magnitude larger than the size of the structures prepared with chlorobenzene solution. According to Nillson et al., [12] who examined the phase behavior in a polyfluorene copolymer:PC60BM mixture, such islands are a result of a high solvent evaporation rate (chloroform, toluene, or xylene) and a strong polymer:fullerene repulsive interaction. Here, our results partly support the first of Nillson's observations. For two highly volatile solvents, toluene and chloroform, we obtained cluster structures only for some mixing ratios. Additionally, in our experiment, the interaction between PCDTBT and PC70BM was relatively low and did not vary within mixtures; nevertheless, both clusters and flat layers were formed. Here, we favor the explanation presented by Michels and Moons. [10] They demonstrated that for systems with a critical point located in the high solvent fraction region, the transition between lamellar and lateral structures can take place for some mixing ratios. Such lamellar structures are visible in depth profiles obtained for films casted from chloroform (blend ratios 1:3 and 1:4) and o-dichlorobenzene (blend ratios 2:1, 1:1, and 1:3). Gradient in the PCDTBT concentration is also visible for sample prepared from chlorobenzene solution, mixing ratio 1:4. In the case of toluene:PCDTBT:PC70BM, we did not observe the critical point, which means that from the beginning, the system is located within an unstable region. For low polymer concentrations, the mobility of the blend components is sufficiently high to allow a drop-like breakup.
Phase separation in o-dichlorobenzene:PCDTBT:PC70BM blends leads to very fine structures that are almost identical. The roughness of the free surface is ca. 0.25 nm for all tested mixture ratios. A similar morphology was also formed for the 2:1 mixture spin-cast from chlorobenzene, but the very fine grain-like structures become more pronounced for higher contents of PC70BM (1:3 and 1:4). The roughness of the o-dichlorobenzene films that were cast varies between 0.25 and 0.45 nm. Long wavelength undulations of the free surface clearly distinguish the layers cast from chloroform solutions from films prepared from the other solvents under study. These structures are also higher than those observed for chlorobenzene and o-dichlorobenzene; here, the roughness reaches a value between 0.8 and 1.3 nm. Our results are compatible with the values presented by Shin et al., [35] who compare the roughness for PCDTBT:PC70BM prepared with chlo-roform, chlorobenzene, and o-dichlorobenzene. They also observed an increase in roughness in the following order: o-DCB to CF to CB.
The ternary phase diagrams for chloroform and chlorobenzene are almost identical, so based only on diagram analysis, the final morphology of the polymer:fullerene layer should be similar. However, there are visible discrepancies in lateral and vertical phase separation. Michels and coworkers showed [10,11] that the dynamics of the phase separation and vitrification, which unfortunately cannot be seen in the phase diagram, are an important factor determining the final morphology. The dynamics are related to the solvent evaporation rate because the faster the evaporation, the more pronounced the lateral structures, [12] which was also qualitatively confirmed in the presented results.

CONCLUSIONS
The first systematic analysis of the morphology of PCDTBT: PC70BM films with different blend ratios cast from different solvents in terms of the ternary phase diagrams was shown. The interaction parameters between PCDTBT and four solvents, namely chloroform, chlorobenzene, o-dichlorobenzene, and toluene, were evaluated from swelling experiments. These values were used to calculate the ternary phase diagrams for PCDTBT:PC70BM:solvent mixtures. The spinodal and binodal lines were obtained for blends dissolved in chloroform, chlorobenzene, and o-dichlorobenzene, whereas for toluene, only the spinodal line was obtained. The lack of a one-phase region in the latter case is due to the strong repulsive interactions between the components of the mixture and the solvent; for this system, the largest variety of structures was observed. Phase separation in o-dichlorobenzene:PCDTBT: PC70BM blends led to very fine structures that were almost independent of the blend composition. The ternary phase diagrams for chloroform and chlorobenzene were almost identical, so based only on diagram analysis, the final morphology of the polymer:fullerene layer should be similar, but there were visible discrepancies in lateral and vertical phase separation. This indicates that the dynamics of the phase separation and vitrification are important factors determining the final morphology.

APPENDIX I: CALCULATION OF THE BIONODAL LINES IN TERNARY SYSTEMS
In general, the Gibbs free energy is a function of the abundances of particles N i , pressure p, and temperature T: N 1 , . . . , N n , p, T) We can assume that in the mixing process, the pressure and temperature are constant, so the chemical potentials can be defined as: The Gibbs free energy is an extensive thermodynamic function, meaning that it scales as G → λG when N i → λN i . It is often easier to work with its intensive version g, defined as: with the numbers x 1 ,..., x n summing up to one: We can express the chemical potentials μ i at constant pressure p and temperature T in terms of g as: x j ∂ j g + ∂ i g or, using a vector notation: E = (1, . . . , 1) Vector X = (x 1 ,…, x n ) satisfies the relation: The binodal is defined as a set of points with equal chemical potentials, that is, a solution that follows the system of equations: Let us now restrict ourselves to ternary mixtures so that . If the system of conditions (A1) does not degenerate (see (A5)), then it imposes 5 constraints on 6 variables, the set of solutions is going to be a curve in a 6-dimensional space, we can treat X 1 and X 2 as functions of some real parameter t. Taking a derivative of (1) with respect to t yields the following set of ordinary differential equations: where and Hess(g)(X) denotes the Hessian of the function g at point X.
Note that H 1 and H 2 are functions of X 1 and X 2 , respectively, so (A2) is a nonlinear system of differential equations for X 1 and X 2 .
If we denote K i = H i · X i ', provided that H i is invertible, the system (A2) simplifies to: In the definition of E i , we used the fact that H i is symmetric, and so is (H i ) −1 . Remarkably, (A3) can be solved explicitly by guessing that the vectors K 1 and K 2 should be both equal to some vector K * that is orthogonal to both E 1 and E 2 . Indeed, if we set K 1 = K 2 = K * = E 1 × E 2 , it is straightforward to see that K 1 and K 2 satisfy (A3). Returning to variables X i ', we obtain a solution to (A2) as follows: provided that the matrices H 1 and H 2 are invertible, that is: Under the condition (A4) and (A2) are equivalent up to a rescaling of K * , we will choose this scaling later on. Suppose now that we have a pair of points, X 1 (0) and X 2 (0), that solve the original binodal condition (A1). By means of (A4), we can find a continuum of solutions X 1 (t) and X 2 (t) by integrating the (nonlinear) system of differential equations (A4), with initial conditions X 1 (0) and X 2 (0).
In the proposed approach, we use the properties of the critical point to define the starting points to solve (A4). We can assume that critical point X C is the point where the two branches of the coexistence curve (X 1 (t) and X 2 (t)) originate. The two curves X 1 (t) and X 2 (t) do not intersect later on, thus X 1 (t) ≠ X 2 (t) for t ≠ 0, and they have to obey the symmetry relation X 1 (t) = X 2 (−t) for all t. Starting with the initial conditions X 1 (0) = X 2 (0) = X C and assuming that the binodal curve is analytic near X C , we can write: Now, plugging this relation back to (A2) and projecting onto R and X C , we obtain the following two conditions: A critical point is a point where a binodal is tangent to a spinodal (a curve where det(H) = 0), defined also as a point X C where: for some vector R. Because such R is tangent to the spinodal at X C and thus to the binodal at the same point, the two binodal curves coming out of it can be expressed as: We will use these two curves as a starting point for solving the system (A4). Note that in (A7), the curves X 1 and X 2 start from the same point X C but go in opposite directions. Moreover, the coefficients in front of R in both X 1 and X 2 are equal up to a point because we expect X 1 and X 2 to be interchangeable after replacing t with −t (i.e., X 1 (−t) = X 2 (t) and vice versa).
We start the integration using the initial data (A7) for some small t, say t = 0.001, and integrate until one of the arms of the binodal reaches an edge of the phase diagram. Such an algorithm produces the two arms of the binodal extending out of the critical point X C . Additionally, before we start the integration, we explicitly use the condition X i · E = 1 to replace one of the coefficients. If we denote the components of X 1 and X 2 as: X 1 = (x 1 , y 1 , z 1 ) , X 2 = (x 2 , y 2 , z 2 ) this means that according to X i · E = x i + y i + z i = 1, we can replace z i with 1 − x i − y i and discard the two equations for z' i . In this way, we reduce the dimension of the differential Eq. (A1) from 6 to 4. Depending on the particular form of g, we can also remove any singular behavior by normalizing K * to one. Another normalization strategy would be to remove any singular terms in X' 1 and X' 2 by multiplying K * by a common denominator of both X' 1 and X' 2 .
The above algorithm is sufficient for qualitatively simple plots containing only a single critical point and no miscibility gap or regions with three coexisting phases.

Electronic Supplementary Information
Electronic supplementary information (ESI) is available free of charge in the online version of this article at http://dx.doi.org/ 10.1007/s10118-020-2424-8.