Pattern formation aspects of electrically charged tri-stable media with implications to bulk heterojunction in organic photovoltaics

Self-organization is a universal property in which spatial and/or temporal order/correlation in systems (or certain symmetry) emerges from interaction of individual subsets at much lower scales [1]. Unlike autocatalytic processes that often give rise to complex temporal dynamics in biological, ecological and chemical systems [2–4], applications related to energy and material science view the emerging self-assembled phenomena as gradient flows [5], i.e., systems that evolve via overall free energy minimization [6]. For the latter, the two fundamental ingredients incorporate the thermodynamic free energy and the potential energy dissipation, here it is convenient to imagine a binary system of two immiscible phases, e.g., oil and water or di-block copolymers, that exhibits asymptotic phase separation. In case the phases are additionally electrically charged, e.g., electrochemical applications, the bulk evolution is driven by competing short-range and long-range interactions and results in rich diversity of periodic lamellar- and hexagonal-like morphologies [7–9]. Moreover, recently it was found that these systems exhibit also several universal properties of dissipative media [10–12], such as formation of coexisting localized states via homoclinic snaking [13], which not only challenges our understanding of electrically charged self-assembly but also has implications to material science and applications thereof.

The ability to design and control over a broad range of meso- and nano-scale morphologies in organic materials makes them attractive to many applications and specifically to renewable energy applications, examples of which include incorporation of ionic liquids in super-capacitors and batteries [14–16], polyelectrolyte membranes used in fuel cells [17–19], and organic photovoltaics (OPV) [20–22]. The latter are particularly intriguing since unlike inorganic solar cells, they are involved with creation and annihilation of an electrically neutral particle —an exciton— that is generated upon illumination [23,24]. Excitons are short-lived charge pairs with a typical diffusion length of roughly 10 nm. Therefore, OPV efficiency critically relies on an intertwined (mosaic) morphology of electron donor and acceptor materials of about exciton diffusion length, enabling the majority of the generated excitons to reach the donor/acceptor interface and dissociate to free charges that are then collected at the electrodes. This fine nano-structure, also known as bulk heterojunction (BHJ), increases the efficiency of the OPV to a level competitive with that of inorganic devices [24–30].

Performance and efficiency of an BHJ OPV depends among many factors also on the integrity of the labyrinthine-type morphology [30–39]. Indeed, even 1 nm changes in the acceptor or donor labyrinth width can lead to a significant decrease in overall efficiency [40,41]. Nevertheless, following the assumption that morphology evolution is much slower than other degradation mechanisms (e.g., oxidation), existing OPV models consider "frozen" morphologies [42–45]. Yet, since exciton dissociation is responsible for electrical charge and, thus, for electrical force, the interrelation between interface properties and charge separation kinetics may become equally important to BHJ evolution, including instability. Furthermore, recent experimental evidences report on identification of an intermediate (third) phase between the donor and the acceptor [46–51], which is often also being referred to as "spaghetti" or "river and streams" phase [52]. The average donor and acceptor concentrations throughout the "spaghetti" phase are about the same. Yet, this intermediate phase is distinct from a randomly distributed amorphous composition, as it posseses correlations of perculative nature and is, therefore, regarded as a third phase of the material. It is believed that properties of this third phase, e.g., composition and width, are strongly related to efficiency of exciton dissociation and also implicitly to the charge transport toward the electrodes [27,52,53]. As such, decoupling between Coulombic forces and the interfacial BHJ evolutions cannot be neglected [42].

Here we extend a continuum (mean-field) framework introduced by Buxton and Clarke [42] by incorporating morphology evolution and accounting for the possible emergence of a stable intermediate phase between the pure donor and acceptor phases. The outcome is a self-consistent model which unifies morphology and electrokinetics by incorporating and coupling three components: multi-interface development employing Onsagers' approach to phase separation, charge transport via modification of Poisson-Nernst-Planck framework, and bimolecular interaction terms for exciton dynamics. Through bifurcation analysis we obtain a parameter plane for coexisting front and periodic solutions involving both two-phases and three-phases, and show that the width of the intermediate phase is in fact independent of the domain size. The analysis is then applied to OPV efficiency in bilayer and ideal stripe morphologies. The approach not only provides a systematic road map to understand BHJ in OPV but also insights into degradation mechanisms related to morphological instabilities and a methodology to tackle three-phase compositions in material science applications.

We start with formulating model equations in which the evolution of donor/acceptor morphology depends on the charge charier concentrations and vice versa. The donor/acceptor morphology is described by an order parameter $u:=\varphi_A-\varphi_D$, where $\varphi_A$ and $\varphi_D$ are volume fractions of acceptor and donor, respectively, so that u varies from $u\equiv u_- = -1$ for pure donor to $u\equiv u_+ = 1$ for pure acceptor. The electrical charges are holes (p) and electrons (n), while in addition there are also neutral excitons $(\chi)$ from which the charges are dissociated. Finally, the system is complemented by Poisson's equation. The modified transport equations are derived from a free energy (see the Supplementary Material Supplementarymaterial.pdf (SM) for details) and following Buxton-Clarke [42] they are supplemented with exciton dissociation/decay kinetics together with electron-hole recombination. Due to the interest in the phenomenological properties of donor/acceptor interfaces, we study the dimensionless version of the model eq. (1); for dimensional forms we refer the reader to the SM,

$$\begin{align} \frac{\partial {u}}{\partial {t}}&=\underbrace{D_u\nabla\cdot{\{{{\nabla u+}{({1-u^2})} {[{{({\beta W''(u)\nabla u-\lambda\nabla^3u})}}]}}\}}}_{\rm phase~separation} +\underbrace{D_u \zeta\nabla\cdot {\{{ {({1-u^2})} {[{{({p+n})}\nabla u+{({1+u})}\nabla p-{({1-u})}\nabla n}]}}\}}}_{\rm donor/acceptor~affinity~to~charges}, \end{align} \tag{ 1a }$$

$$\begin{align} {\frac{\partial {\chi}}{\partial {t}}}&=\underbrace{\nabla^2 \chi}_{\rm diffusion}-\underbrace{\frac{1-u^2}{\tau}\chi}_{\rm dissociation}-\underbrace{\chi}_{\rm decay}+\underbrace{G}_{\rm generation}, \end{align} \tag{ 1b }$$

$$\begin{align} {\frac{\partial {p}}{\partial {t}}}&=D_p\nabla\cdot{[{\underbrace{p\nabla\phi+\nabla p}_{\rm drift-diffusion}+\underbrace{\zeta p(1+u)\nabla u}_{\rm charge~affinity}}]}+\frac{1-u^2}{\tau}\chi-\underbrace{\gamma\,np}_{\rm recombination}, \end{align} \tag{ 1c }$$

$$\begin{align} {\frac{\partial {n}}{\partial {t}}}&=D_n\nabla\cdot{[{-n\nabla\phi+\nabla n-\zeta n(1-u)\nabla u}]}+\frac{1-u^2}{\tau}\chi-\gamma\,np, \end{align} \tag{ 1d }$$

$$\begin{align} &\hspace{-1pc}\nabla \cdot {[{\epsilon \nabla \phi}]}=n-p. \end{align} \tag{ 1e }$$

The first term in (1a) is a typical Cahn-Hilliard approach to phase separation (in the absence of charges) but here we use a fourth-order potential, $W(u)=(1-u^2)(u^2+\xi)/2$, which together with entropy (see the SM) allows the transition from double- to triple-well in the free-energy landscape, i.e., at $\xi<\xi_0=(1+\beta)/\beta$ for the latter. We refer to the additional third phase as an intermediate phase, $u\equiv u_0$ (here $u_0=0$). Keeping fidelity to BHJ OPV, electrical charges originate from exciton dissociation, whose rate is maximal at the donor/acceptor interface, and specifically at $u=u_0$, while being negligible at the pure phases, $u=\pm1$ [42]. Thus, the dissociation term is proportional to $1-u^2$, as seen in (1b). Due to presence of electrical charges, the length scale of the intermediate phase is also affected by the second term in (1a) which represents coupling to electrons and holes, similarly to the interaction in the context of charged polymers [54]. The other eqs. (1c) and (1d) comprise standard terms with additional contribution to the transport of charge preference to respective phases, i.e., charge affinity. Model (1) not only captures all the qualitative features of charge generation-recombination and transport that have been suggested by Buxton and Clarke [42] but also allows coupling between Coulombic forces to morphology evolution. Notably the introduction of the intermediate phase does not require using Onsager's phenomenological theory for dissociation and the Poole-Frenkel mobility, that have been used to evaluate charge dynamics over a frozen double-well media [42]. This distinct approach provides a framework that is not limited to specific functional forms, the formulation is flexible and allows physical refinement by demand. The qualitative features of one- and two-space dimensional solutions, however, will be shown as generic since they rely on domain walls (or what is mathematically being referred to as heteroclinic connections), so that the quantitative regimes can change while the qualitative properties persist.

To understand the pattern formation mechanisms and implications of (1), we seek first solutions in one spatial dimension (1D) and then relate them to a more realistic OPV-type configuration in 2D. We choose to focus on two control parameters: the first is G > 0 which is related to exciton generation due to illumination and hence to charge density, and the second is ξ that corresponds to the depth of the intermediate well.

Equation (1) have a multiplicity of uniform solutions due to conservation of the order parameter u, see (1a), which are given by

$$\begin{equation} \left(\begin{array}{c} u_0\\\chi_0\\p_0\\n_0\\\phi_0 \end{array}\right)=\left(\begin{array}{c} u\\[5pt] \displaystyle\frac{\tau G}{\tau+1-u^2}\\[10pt] \sqrt{\displaystyle\frac{G}{\gamma}\displaystyle\frac{1-u^2}{\tau+1-u^2}}\\[10pt] \sqrt{\displaystyle\frac{G}{\gamma}\displaystyle\frac{1-u^2}{\tau+1-u^2}}\\[10pt] 0 \end{array}\right), \end{equation} \tag{ 2 }$$

and determined by selection of u. As for the typical bi-stable system (double-well potential, $\xi>\xi_0$) also here the pure phases $u= |u_\pm| \lesssim 1$ are linearly stable while u₀ is linearly unstable for any G > 0. Stability properties of u₀ for the tri-stable case, i.e., for $\xi<\xi_0$, are different as will be shown next.

As expected linear stability of pure phases $u=u_\pm$ persists everywhere and thus the existence of front solutions (under Neumann boundary conditions), i.e., solutions for which $u \to u_\pm$ as $x\to \pm \infty$. On the other hand, the linear stability analysis about u₀ of the form [10,11,13]

$$\begin{equation} \left(\begin{array}{c} u\\\chi\\p\\n \end{array}\right)-\left(\begin{array}{c} 0\\\chi_0\\p_0\\n_0 \end{array}\right)\propto \text{e}^{\sigma t+\text{i}kx}+{\rm complex~conjugated}, \end{equation} \tag{ 3 }$$

with substitution of the linearized Poisson's equation into other fields, shows that one dispersion relation corresponds to a finite wave number instability, namely at $G=G_c$ the growth rate σ of the perturbations is real and satisfies $\sigma(0)=\sigma(k_c)=0$ and otherwise $\sigma(k)<0$, for all k. This instability is, however, coupled to the marginal Goldstone mode k = 0 due to the conservation of the order parameter, i.e., $\sigma(0)=0$ always persists. The right inset in fig. 1(a) shows a typical dispersion relation for G > 0 and $\xi<\xi_0$.

Zoom In Zoom Out Reset image size

While the dispersion relations are cumbersome and obtained numerically, the critical values at the onset can be derived:

$$\begin{align} G_c&=\frac{\epsilon^2\gamma(\tau+1)(\beta\xi-\beta-1)^2}{{4(\lambda-\epsilon\zeta+\epsilon\zeta^2-2\zeta\sqrt{\epsilon\lambda})^2}}, \end{align} \tag{ 4 }$$

$$\begin{align} k^2_c&=\frac{(1-\zeta\sqrt{\epsilon/\lambda})(\beta\xi-\beta-1)}{{\lambda-\epsilon\zeta+\epsilon\zeta^2-2\zeta\sqrt{\epsilon\lambda}}}. \end{align} \tag{ 5 }$$

Consequently, the onset properties within this model are found to be independent of diffusion coefficients. On the other hand, while the onset G_c depends linearly on the recombination rate $(\gamma)$ and the dissociation rate $(\tau)$ there is no dependence on the periodicity, i.e., change in k_c. The latter indicates that excitons play a negligible role in periodicity of the bifurcating periodic solutions, u_p, which are determined mostly by material properties and interaction between the material and charges, i.e., affinity.

Periodic solutions (u_p) bifurcate sub-critically from $u=u_0$ at the onset $G=G_c$, that is toward the linearly stable portion of u₀, and thus are unstable, as demonstrated in fig. 1(b). They gain stability via a saddle-node bifurcation about $G=G_{SN}$ and fold back to large values of G. The profiles about the intermediate phase of this stable periodic solutions coincide with those of the coexisting front solutions (u_f) that are superimposed in fig. 1(b). This similarity implies a generic property that we recover using spatial dynamics analysis, i.e., omitting temporal derivatives, rewriting (1) as first-order ordinary differential equations, and looking at steady-state solutions and envision space as time in a similar fashion as for the Ohta-Kawasaki analysis used for charged polymers [13]. Linearization about u₀ reveals twelve eigenvalues $(\mu)$ of which quartic multiplicity is always at zero (Re μ = Im μ = 0), four are always real, and the remaining four play a crucial role determining the signature of the intermediate phase in profiles that are shown in fig. 1(b). At the finite wave number instability onset $(G=G_c)$, these four eigenvalues lie in a double multiplicity and are purely imaginary (Re μ = 0) which is a signature of a reversible Hopf bifurcation (in space) [13,57,58], as shown in fig. 2. For G < G_c these eigenvalues split and for G_T < G < G_c persist as complex conjugated, implying the connection to the unstable manifold of unstable limit cycles or a periodic solutions in space that bifurcate from $G=G_c$, as shown in the top panel of fig. 2. At $G=G_T$, where

$$\begin{equation} G_T=\frac{\epsilon^2\gamma(\tau + 1)(\beta\xi-\beta-1)^2}{4(\lambda-\epsilon\zeta+\epsilon\zeta^2+2\zeta\sqrt{\epsilon\lambda})^2}, \end{equation} \tag{ 6 }$$

the eigenvalues collide on the real axis and since all eigenvalues are real the impact of the limit cycles vanishes, namely for G < G_T the attractor about u₀ is hyperbolic in space and thus the orbit spends more "time" near u₀ making the inflection in the profile visible, as is also shown in the top panel of fig. 2. The analytical result for G_T (dashed line in fig. 1(a)) agrees well with the numerical computations (dotted line). However, the situation is not necessarily generic and under different parameter setting $(G_T \ll G_c)$, the spatial structure about the u = 0 state can be expected to show also decayed oscillations in space [59], yet, this is beyond the scope here and will be addressed elsewhere. Notably, although mathematically it is obvious that G_T < G_c due to temporal linear stability, the result is also consistent physically: condition G_T > G_c implies that reference energies are lower than the thermal energy $(\zeta<k_BT)$ which is an unrealistic inequality, see the SM for details and definitions.

Zoom In Zoom Out Reset image size

To relate the 1D analysis to the role of u₀ in certain aspects of OPV performance, we consider a two-dimensional domain $[x\times y]\in [L_x \times L_y]$ and also allow the charge collection from the electrodes that are placed along the y boundaries (see fig. 3(a)). The overall boundary conditions for eqs. (1) include periodic boundary conditions in the x-direction, charge flux (p and n) and fixed potential in the y-direction [60]:

$$\begin{equation*} (J^u_y,J^\chi_y,J^p_y,J^n_y,\phi)\vert_{y=0}=\biggl(0,0,-D_p p\frac{\partial \phi}{\partial y},0,\frac{V}{2}\biggr) \end{equation*}$$

and

$$\begin{equation*} \left(J^u_y,J^\chi_y,J^p_y,J^n_y,\phi \right)\big\vert_{y=L_y}=\left(0,0,0,-D_n n\frac{\partial \phi}{\partial y},-\frac{V}{2}\right), \end{equation*}$$

where J_y for each field is given in the SM, V is the applied voltage between the electrodes which here we take to be under short circuit conditions, i.e., V = 0.

Zoom In Zoom Out Reset image size

First, we examine the effect of the applied voltage on a bilayer geometry, i.e., solution u_f with respect to the y-direction (fig. 3(a)). Change in the boundary condition for u_f results in the current of charges (finite-size effect) which, in turn, increases Λ, although approaching u_f that is obtained in one space dimension (i.e., without charge collection) when $L_y \to \infty$, as demonstrated in fig. 3(b). The increase in Λ is profound for $G \ll G_T$ and becomes more moderate as $G \to G_T$. In fact, the increase in Λ becomes crucial for ideal (stripe) morphology, as shown in fig. 3(c); stripe morphology is obtained by extending the periodic solutions (u_p) in y and to keep fidelity to OPV BHJ, we introduce an additional layer of the phases near the electrodes.

The increase in Λ results in an instability of stripes so that stripes involving u₀ persist only near G_T, where $\Lambda\ll\ell_c=2\pi/k_c$; here from about $G\approx3$ and up to G_T, while for G > G_T the intermediate phase, u₀, vanishes $({\Lambda=0})$ as in the one-dimensional case. In particular, in this regime the quantitative analogy between the bilayer and stripes with respect to Λ still persists, thus the width of the intermediate phase, u₀, in a 2D device under applied voltage can be efficiently deduced from the analysis in 1D, either from u_f or u_p.

As expected, the trend of the current density, $j^n_y=L_x^{-1}\int J^n_y{\rm d}x$, at short circuit conditions increases with G (see inset in fig. 3(b)) while the ratio $j_y^n/\bar{n}$, where $\bar{n}=L_x^{-1}L_y^{-1}\int n(x,y) {\rm d}x {\rm d}y$, that indicates outflow vs. creation (due to G) of charge is essentially constant. In fact, this region about G_T of constant normalized current indicates the stability of periodic solutions (ideal BHJ), a direct consequence of the tri-stable system. Specifically, outside this region, i.e., at higher G values, periodic stripes are not influenced by the attractor near u₀, and become unstable to zig-zag [13] leading to deformed morphology, as shown in fig. 4. Notably, a detailed analysis of stripe instabilities is beyond the scope of this letter and will be discussed in detail elsewhere, e.g., a possible transition to isolated domains due to varicose instability that decreases the cell's efficiency [41].

Zoom In Zoom Out Reset image size

We have formulated a distinct mean-field model that enables to capture binary (bi-stable) and ternary (tri-stable) bulk heterojuction (BHJ). We showed that although the width of the intermediate phase $(\Lambda)$ in BHJ depends on illumination strength (G) and applied voltage (V), its trend can be efficiently deduced from a relatively simply analysis in one space dimension as long as $L_y\gg \Lambda$, which is consistent with device-wise applications. Moreover, we conjecture that the existence of the intermediate phase impacts stability of BHJ striped morphology —a criterion that is valuable to design efficient and stable BHJ-based organic photovoltaic devices. In a broader context of material science, our framework with appropriate modifications of the free-energy contractions, opens new vistas to a wider range of electrically charged amphiphilic-type mixtures to tackle their stability and evolution, examples of which include emulsions, polyelectrolytes, perovskite films, and ionic liquids.

We thank Prof. Iris Visoly-Fisher and Prof. Rafi Shikler (Ben-Gurion University of the Negev) for stimulating and helpful discussions on OPV and Prof. Hannes Uecker (Universität Oldenburg) for help in the implementation of the pde2path package. The research was done in the framework of the Grand Technion Energy Program (GTEP) and of the BGU Energy Initiative Program, and supported by the Adelis Foundation for renewable energy research.