Bending and pinching of three-phase stripes: From secondary instabilities to morphological deformations in organic photovoltaics
BBending and pinching of three-phase stripes: From secondaryinstabilities to morphological deformations in organic photovoltaics
Alon Z. Shapira, Nir Gavish, Hannes Uecker, and Arik Yochelis
4, 5, ∗ Swiss Institute for Dryland Environmental and Energy Research,Blaustein Institutes for Desert Research,Ben-Gurion University of the Negev, Sede Boqer Campus,Midreshet Ben-Gurion 8499000, Israel Department of Mathematics, Technion - IIT, Haifa, 3200003, Israel Institute for Mathematics, Carl von Ossietzky University of Oldenburg,P.F 2503, 26111 Oldenburg, Germany Department of Solar Energy and Environmental Physics,Blaustein Institutes for Desert Research,Ben-Gurion University of the Negev, Sede Boqer Campus,Midreshet Ben-Gurion 8499000, Israel Department of Physics, Ben-Gurion University of the Negev, Beer Sheva 8410501, Israel (Dated: November 17, 2020) a r X i v : . [ n li n . PS ] N ov bstract Optimizing the properties of the mosaic morphology of bulk heterojunction (BHJ) organic pho-tovoltaics (OPV) is not only challenging technologically but also intriguing from the mechanisticpoint of view. Among the recent breakthroughs is the identification and utilization of a three-phase(donor/mixed/acceptor) BHJ, where the (intermediate) mixed-phase can inhibit morphological changes,such as phase separation. Using a mean-field approach, we reveal and distinguish, between generic mech-anisms that alter through transverse instabilities the evolution of stripes: the bending (zigzag mode)and the pinching (cross-roll mode) of the donor/acceptor domains. The results are summarized in a pa-rameter plane spanned by the mixing energy and illumination, and show that donor-acceptor mixtureswith higher mixing energy are more likely to develop pinching under charge-flux boundary conditions.The latter is notorious as it leads to the formation of disconnected domains and hence to loss of chargeflux. We believe that these results provide a qualitative road-map for BHJ optimization, using mixed-phase composition and therefore, an essential step toward long-lasting OPV. More broadly, the resultsare also of relevance to study the coexistence of multiple-phase domains in material science, such as inion-intercalated rechargeable batteries. ∗ [email protected] . INTRODUCTION Organic photovoltaics (OPV) are being subjected to intensive research over the past twodecades not only due to their potential advantages as portable and/or lightweight technologicaldevices but also for their intriguing physicochemical mechanisms of operation [1–4]. At the heartof the OPV is the nano-scale mosaic active layer of electron donor (D) and electron acceptor (A)materials, the so-called bulk heterojunction (BHJ) [5–11]. This subtle morphology is essential forefficient dissociation of excitons at the D/A interfaces to electrons and holes and for transportof the latter toward the collectors [12–14]. The short lifetime of the excitons is translated toa spatial length scale, also known as the diffusion length , which respectively sets about tens ofnanometer bi-continuous (ideally comb-like) morphology [11, 15–24].Recent evidences however, indicate that in some compositions [4, 25–33] a third phase, whichis being referred to as a mixed-phase (MP), may additionally become stable along with the pureD/A phases [34, 35]. This MP has a molecular percolating structure about a 1:1 ratio between thedonor and the acceptor molecules [35] and thus is distinct from a random distribution although inboth cases the averaged quantity is identical. As such, the MP can be thought of as an effectiveenergetic barrier (as being an intermediate metastable state) between the energetically favorableD and A phases [35, 36] while keeping the exciton dissociation properties intact. Recent studiesindicate that MP plays a role in the evolution of BHJ, ranging from the width and form of theD/A interface [25, 36, 37] to an inhibitor of the phase separation process [35].(a) (b)
Figure 1. (a) Schematic illustrations of (a) the bending (ZZ) and (b) the cross-roll (CR) instabilities.Horizontal arrows indicate the bending direction while vertical arrows show the directions of the oppositecharge fluxes to the electrodes.
II. DETERMINING THE DONOR-ACCEPTOR RATIO
In the dark, the free energy comprises the entropy and the mixing energy for the materialorder parameter [36], u := ϕ A − ϕ D ∈ [ − , , where ϕ A , ϕ D are the respective fractions of theA/D phases. Its dimensionless form reads E M ( u ) = (cid:90) Ω − u − u (cid:124) (cid:123)(cid:122) (cid:125) donor’s entropy + 1 + u u (cid:124) (cid:123)(cid:122) (cid:125) acceptor’s entropy + β − u )( u + ξ ) + λ |∇ u | − e (cid:124) (cid:123)(cid:122) (cid:125) mixing energy d x , (1)where Ω is the domain, which we take to be a rectangle Ω = (0 , l x ) × (0 , l y ) . Further, e is areference energy density for which the minimum of E M ( u ) is zero, β determines the ratio betweenmixing energy and entropy, and ξ determines the depth of the intermediate well such that small ξ corresponds to a lower mixing energy (see Fig. 2(a)), and λ is the penalty for creation ofmultiple interfaces and associated with the width of the interface. Due to entropy, the minimumenergy of donor-rich ( u := u − ) and acceptor-rich ( u := u + ) phases are shifted from u = ± to slightly lower values in | u | , while the mixed-phase always sits at u := u = 0 , as shown inFig. 2(a). Due to conservation of the order parameter, however, there are many other uniformsolutions u = u ∗ and these solutions are related to the D:A ratio of non-uniform solutions, e.g.,D-A interfaces. The connection between the u ∗ and the D:A ratio is made through averaging of4 in one space dimension (1D), (cid:104) u (cid:105) := l − x (cid:90) l x u d x. (2)For the symmetric case (cid:104) u (cid:105) = 0 , the amount of donor and acceptor is identical so that theinterface is located at x = l x / , and for the asymmetric case, where |(cid:104) u (cid:105)| > , this location isshifted; note that for the uniform states (cid:104) u (cid:105) = u ∗ . Thus, for non-uniform solutions that are ofinterest here, it is required to identify the allowed range of (cid:104) u (cid:105) and we do it by looking at thestability of u ∗ .The evolution equation [36] (in the dark and with mobility D u (1 − u ) ) reads as ∂u∂t = D u ∂ u∂x + D u ∂∂x (cid:26)(cid:0) − u (cid:1) (cid:20) β (1 − u − ξ ) ∂u∂x − λ ∂ u∂x (cid:21)(cid:27) , (3)where D u is the diffusion coefficient. Linear stability analysis (performed on an infinite domain)about uniform states u = u ∗ corresponds to u − u ∗ ∝ e αt + ikx + c.c. , (4)where c.c. is the complex conjugate and α is the perturbation growth rate of wavenumber k andgiven by α ( k ) = − D u k (cid:8) − u ∗ ) (cid:2) β (1 − u ∗ − ξ ) + λk (cid:3)(cid:9) . (5)The instability of u = u ∗ is of a typical long-wavenumber type [40] and the regime of unstablesteady state solutions u min ∗ < | u ∗ | < u max ∗ is obtained by taking the limit α ( k ) → as k → (seeFig. 2(b)), where u min ∗ = 12 √ (cid:113) − ξ − (cid:112) (5 + ξ ) − /β, u max ∗ = 12 √ (cid:113) − ξ + (cid:112) (5 + ξ ) − /β. These results imply that any choice within the u min ∗ < | u ∗ | < u max ∗ range will result in a coarseningdynamics that minimizes the energetic penalty of interfaces, following which phase separation isachieved (not shown here). In the inset of Fig. 2(b), we show that indeed the interface solutionsare bi-asymptotic to u ± and that the location of the interface shifts according to (cid:104) u (cid:105) that is setby u ∗ . Thus, in the analysis that follows, we will focus on the regime < | u ∗ | < u min ∗ , underillumination. 5a) (b) Figure 2. (a) Free energy functional (1) for uniform u , with ξ = 2 . , e = − . (blue, genuine triplewell) and ξ = 2 . , e = 0 . (red, closer to double well), where ξ corresponds to the mixing energy.(b) Stability and instability intervals for uniform solutions determined by the critical values u min ∗ ≈ . and u max ∗ ≈ . . The inset demonstrates three-phase interface solutions, as computed from (3), for (cid:104) u (cid:105) = 0 and (cid:104) u (cid:105) = 0 . (as marked by ‘ • ’ in the main figure, respectively); the location of the interfacewith (cid:104) u (cid:105) = 0 is at x = 0 . The horizontal solid and dashed lines indicate the value u ∗ = (cid:104) u (cid:105) = 0 . and u ∗ = (cid:104) u (cid:105) = 0 , respectively. Parameters: D u = 1 , β = 0 . , ξ = 2 . , λ = 0 . , k = 10 − . III. EXISTENCE OF STRIPES UNDER ILLUMINATION AND THE EFFECT OF D-A ASYMMETRY
Illumination leads to the creation of excitons which dissociate at the D/A interface and thus,drive the BHJ out of equilibrium. The (dimensionless) total free energy under illumination takesthe form [36] E = E M + (cid:90) Ω χ ln χ + p ln p + n ln n (cid:124) (cid:123)(cid:122) (cid:125) charges entropy + φ ( p − n ) − (cid:15) |∇ φ | (cid:124) (cid:123)(cid:122) (cid:125) electrostatic energy + ζ (cid:2) p (1 + u ) + n (1 − u ) (cid:3)(cid:124) (cid:123)(cid:122) (cid:125) charge affinity d x , where the corresponding equations of motion also incorporate the generation/recombinationfollowing Buxton and Clarke [41], and fluxes of electrical charges coupled to morphological6volution of the BHJ order parameter [36]: ∂u∂t = D u ∇ u + D u ∇ · (cid:8)(cid:0) − u (cid:1) (cid:2) β (1 − u − ξ ) ∇ u − λ ∇ u (cid:3)(cid:9)(cid:124) (cid:123)(cid:122) (cid:125) phase separation + D u ζ ∇ · (cid:8)(cid:0) − u (cid:1) [( p + n ) ∇ u + (1 + u ) ∇ p − (1 − u ) ∇ n ] (cid:9)(cid:124) (cid:123)(cid:122) (cid:125) donor / acceptor affinity to charges , (6a) ∂χ∂t = ∇ χ (cid:124)(cid:123)(cid:122)(cid:125) diffusion − τ − (cid:0) − u (cid:1) χ (cid:124) (cid:123)(cid:122) (cid:125) dissociation − χ (cid:124)(cid:123)(cid:122)(cid:125) decay + G (cid:124)(cid:123)(cid:122)(cid:125) generation , (6b) ∂p∂t = D p ∇ · [ p ∇ φ + ∇ p (cid:124) (cid:123)(cid:122) (cid:125) drift − diffusion + ζp (1 + u ) ∇ u (cid:124) (cid:123)(cid:122) (cid:125) charge affinity ] + τ − (cid:0) − u (cid:1) χ − γ np (cid:124)(cid:123)(cid:122)(cid:125) recombination , (6c) ∂n∂t = D n ∇ · [ − n ∇ φ + ∇ n − ζn (1 − u ) ∇ u ] + τ − (cid:0) − u (cid:1) χ − γ np, (6d) ∇ · [ (cid:15) ∇ φ ] + p − n. (6e)Here the fields χ, p, n stand for excitons, holes and electrons, respectively, φ is the electricpotential, D p , D n are the respective diffusion constants, ζ is the interaction energy between elec-tron/holes and donor/acceptor, τ is the excitons dissociation time, G is the excitons generationrate, γ is the electron-hole recombination rate, and (cid:15) is the permittivity. For details we refer thereader to Shapira et al. [36].Uniform solutions of system (6) are given by U ∗ = ( u ∗ , χ ∗ , p ∗ , n ∗ , , where χ ∗ = τ G/ ( τ + 1 − u ∗ ) and p ∗ = n ∗ = (cid:112) G (1 − u ∗ ) / ( γ ( τ + 1 − u ∗ )) . Linear analysis in 1D, by replacing u ( x ) with U ( x ) , shows that in range ≤ | u ∗ | < u min ∗ , the uniform solution U ∗ goes through a subcriticalfinite wavenumber instability at G = G c , giving rise to periodic solutions U (cid:96) ( x ) with wavenumber k c , that corresponds to the spatial wavelength (cid:96) c = 2 π/k c , where for u ∗ = 0 we get G c = (cid:15) γ ( τ + 1)( βξ − β − λ − (cid:15)ζ + (cid:15)ζ − ζ √ (cid:15)λ ) , and k c = (1 − ζ (cid:112) (cid:15)/λ )( βξ − β − λ − (cid:15)ζ + (cid:15)ζ − ζ √ (cid:15)λ , while for | u ∗ | > the critical values are computed numerically. The periodic solutions U (cid:96) ( x ) =( u (cid:96) ( x ) , χ (cid:96) ( x ) , p (cid:96) ( x ) , n (cid:96) ( x ) , φ (cid:96) ( x )) bifurcate toward the stable portion of U ∗ , that is in direction G < G c , and thus, are initially unstable. Then they grow in amplitude and stabilize after the7a) (b) Figure 3. (a) Bifurcation diagram showing the branches of symmetric (cyan) and asymmetric (blue)periodic solutions in 1D, u (cid:96) , computed for system (6); solid/dashed lines indicate stable/unstable solu-tions. For these computations we employed the numerical continuation package pde2path [42, 43], withperiodic BC. The periodic solutions bifurcate from the uniform solutions u ∗ = 0 and u ∗ = 0 . at therespective values of G = G c . (b) Asymmetric stripe pattern plotted at G = 2 after extending u (cid:96) in the y direction. The solid horizontal line indicates (cid:104) u (cid:105) = 0 . . The arrows at u = ± . indicate the width ofthe acceptor (A) and donor (D) phases, with A being wider. The domain size is Ω = [0 , (cid:96) c ] × [0 , and the critical wavelength of u (cid:96) is (cid:96) c ≈ . . Parameters: β = 0 . , ξ = 2 . , λ = 0 . , ζ = 4 , (cid:15) = 0 . , τ = γ = 100 , D p = D n = 3 , D u = 1 . saddle node bifurcation that is located close to G = 0 , and continue to be stable as G increases,as shown in Fig. 3(a)].Notably, conservation of the order parameter u forces also the periodic solutions u (cid:96) to keepthe average value that is initially set by u ∗ . Namely, periodic solutions (which can be extendedin y direction to form stripes) that bifurcate from u ∗ = 0 correspond to symmetric stripes (i.e.,identical width of the donor and the acceptor domains) while periodic solutions that bifurcatefrom u ∗ = 0 . , for example, are asymmetric , in which acceptor domains are wider; the latteris demonstrated in 3(b). Next, we calculate the stability properties of stripes in the transversedirection. IV. TRANSVERSE INSTABILITY OF THREE-PHASE DONOR/MIXED/ACCEPTORSTRIPES
For the linear transverse instability of stripes to zigzag (ZZ) that corresponds to bendingand to cross-roll (CR) that causes pinching (see Fig. 1), we employ a general space dependenteigenvalue method [44–50] that has been used for example in the context of convection rolls,8tripes in reaction-diffusion media, and thin fluid films. However, due to application to OPV,our interest here is to reveal the impact of physical boundary conditions on the instability ofstripes, i.e., on non-periodic domains in y directions. A. Linear analysis on unbounded domains
We start however, by performing a general analysis of stripes on non-physical infinite (periodicin y direction) domains [39] U ( t, x, y ) − U (cid:96) ( x ) ∝ ˜ U ( x ) e ηt + ik y y + c.c. , (7)where η is the growth rate of the wavenumber, k y , in the transverse direction to U (cid:96) ( x ) , and ˜ U ( x ) = (˜ u ( x ) , ˜ χ ( x ) , ˜ p ( x ) , ˜ n ( x ) , ˜ φ ( x )) is always the periodic eigenfunction. This formulationintroduces a generalized eigenvalue system η M ˜ U = L ˜ U . (8)In (8) M is a singular projection matrix [38] M = , (9) L is a linear operator L [ U (cid:96) ; k y ] = D u L , D u L , D u L , χu/τ L , L , (1 − u ) /τ L , − γp L , L , (1 − u ) /τ − γn L , L , − (cid:15) ( ˆ ∂ x − k y ) , (10)9here L , = λ (cid:2) − u (cid:3) ∂ x +2 λuu x ∂ x + L (0)1 , + L (1)1 , ∂ x + L (2)1 , ∂ x , L (0)1 , = (cid:8) λ (cid:2) uu xxxx +2 u x u xxx − k y (1 − u ) (cid:3) + β (cid:2) u xx (24 u − u )+ u x (72 u − − k y (6 u − u +1)+ ξ (2 u x +2 uu xx + k y (1 − u )) (cid:3) + ζ (cid:2) p xx (1 − u − u )+ n xx (1+2 u − u ) − uu xx + u x )( p + n )+ u x ((1 − u ) n x − (1+4 u ) p x ) − k ( p + n )(1 − u ) (cid:3) − k y (cid:9) , L (1)1 , = (cid:8) λu (cid:2) u xxx − k y u x (cid:3) +4 βuu x (cid:2) u − ξ (cid:3) − ζ (cid:2) uu x ( p + n )+2 p x (2 u − u )+2 n x (2 u − − u ) (cid:3)(cid:9) , L (2)1 , = (cid:8) β (cid:2) u − u +1 − ξ (1 − u ) (cid:3) + ζ ( p + n )(1 − u )+2 λk y (1 − u ) (cid:9) , L , = ζ (cid:2) (1+ u )(1 − u ) (cid:3) ∂ x − ζu x (cid:2) u + u − (cid:3) ∂ x − ζ (cid:8) uu x − u xx (1 − u )+ k y (1+ u )(1 − u ) (cid:9) , L , = ζ (cid:2) ( u − − u ) (cid:3) ∂ x − ζu x (cid:2) u − u − (cid:3) ∂ x − ζ (cid:8) uu x − u xx (1 − u )+ k y ( u − − u ) (cid:9) , L , = ∂ x +1 − k y +(1 − u ) /τ , L , = D p ζp (1+ u ) ∂ x + D p ζ (cid:2) p x (1+ u )+2 pu x (cid:3) ∂ x + (cid:8) D p ζ (cid:2) pu xx + p x u x − k y p (1+ u ) (cid:3) − χu/τ (cid:9) , L , = D p ∂ x + D p (cid:2) φ x + ζu x (1+ u ) (cid:3) ∂ x + (cid:8) D p (cid:2) φ xx − k y + ζ [ u xx (1+ u )+ u x ] (cid:3) − γn (cid:9) , L , = D p (cid:8) p ˆ ∂ x + p x ˆ ∂ x − k y p (cid:9) , L , = D n ζn ( u − ∂ x + D n ζ (cid:2) n x ( u − nu x (cid:3) ∂ x + (cid:8) D n ζ (cid:2) nu xx + n x u x − k y n ( u − (cid:3) − χu/τ (cid:9) , L , = D n ∂ x + D n (cid:2) − φ x + ζu x ( u − (cid:3) ∂ x + (cid:8) D n (cid:2) − φ xx − k y + ζ [ u xx ( u − u x ] (cid:3) − γp (cid:9) , L , = − D n (cid:8) n ˆ ∂ x + n x ˆ ∂ x − k y n (cid:9) . In addition, we employ in (10) spatial operators ∂ x := G , ∂ x := D with periodic boundaryconditions [51] G ≈ x − − − . . . − − , D ≈ x − − − . . . − − , where empty entries are zeros and ∆ x is the spatial distance between two points on the uniformgrid and respectively, the operators ˆ ∂ x := ˆ G , ˆ ∂ x := ˆ D with two-sided homogeneous Dirichlet10a) (b)(c) (d) Figure 4. (a) Dispersion relations showing the growth rates of the unstable ZZ ( η ZZ , green) and the stableCR ( η CR , orange) modes at ξ = 2 . ; at k max y , η ZZ is maximal. (b) Same as (a) but for ξ = 2 . wherealso the CR mode is unstable; here, k max y marks the maximum of η CR . (c) Normalized eigenfunctions ofthe ZZ instability, ˜ u ZZ ( x ) centered around u = 0 and rescaled to one half the amplitude of the periodicsolution, u (cid:96) (gray line) for ξ = 2 . . (d) Same as (c) but for the CR eigenfunction ˜ u CR ( x ) (orange) at ξ = 2 . . Parameters: G = 8 and other parameters as in Fig. 3. boundary conditions to eliminate potential jumps ˆ G ≈ x − − . . . − − , ˆ D ≈ x − − − . . . − − , where, for higher-order derivatives we used the identities ∂ x = GD and ∂ x = D .In Fig. 4, we show two numerical realizations that produce the instabilities schematicallydepicted in Fig. 1, for different MP well depths ξ = 2 . (in (a)) and ξ = 2 . (in (b)) whilekeeping illumination fixed, G = 8 . The dispersion relations ( η ZZ and η CR ) indicate that whilein both cases the ZZ (odd symmetry) mode is unstable for G > G ZZ , the CR (even symmetry)11 igure 5. Parameter plane ( G, ξ ) , corresponding to illumination ( G ) and mixing energy ( ξ ), showingthe instability regions of stripes. Solid/dashed lines denote the onsets G ZZ and G CR for the symmet-ric/asymmetric D-A ratio ( (cid:104) u (cid:105) = 0 / (cid:104) u (cid:105) = 0 . ) as obtained from (8) while ‘ • ’ and ‘ (cid:72) ’ correspond toonsets obtained by numerical integration of (6) with periodic boundary conditions in the x directionand charge flux in the y direction [36] (see text for details). The insets show snapshots of u obtainedby numerical integration of (6): left inset is asymptotic solution at ( G, ξ ) = (2 , . , middle inset for ( G, ξ ) = (8 , . and t = 2800 as indicated by the bottom diamond, and right inset for ( G, ξ ) = (8 , . and t = 1200 as indicated by the top diamond. The green and orange envelope lines in the middle andthe right insets represent the ZZ and the CR modes, respectively, as obtained by eigenvalue analysis (8):the modes are parameterized as x = (4 ∓ / (cid:96) c − ε cos( k max y y ) with ε = 0 . and k max y = 0 . as in-dicated in Figure 4(a) and x = (4 ∓ / (cid:96) c ± ε cos( k max y y ) with ε = 0 . and k max y = 0 . as indicatedin Figure 4(b), respectively. The simulations were performed on a domain Ω = [0 , (cid:96) c ] × [0 , : forZZ with (cid:96) c ≈ . in x at ξ = 2 . and for CR with (cid:96) c ≈ . at ξ = 2 . . Other parameters are as inFigure 3. mode becomes unstable (with η CR ( k max y ) > ) only above G = G CR ; in (c,d) we also show thecorresponding ˜ u component of the eigenfunctions for the ZZ or the CR modes at k max y .We generalize the results in a parameter plane spanned by ( G, ξ ) (see Fig. 5) and show that asimilar trend persists (dashed lines) also for the asymmetric donor-acceptor ratio. The instability12 nsets are defined such that the maximal growth rate becomes positive, i.e. when η ( k max y ) > ,at G ZZ and G CR , respectively. This implies degeneracy for G > G CR , the onset of CR mode, aregion in which a competition between bending and pinching of stripes should be expected (eventhough the ZZ mode has a larger growth rate, η ZZ > η CR ). Next, we show that this degeneracyis destroyed once we allow passage of current through boundaries in the y direction, i.e., physicalboundary conditions. B. Realization of instability modes in the presence of charge outflux
Although in the above analysis we used non-physical boundary conditions as we did not allowcharge flux through the boundaries in y direction, the results provide a good guiding for realisticcharge-flux boundary conditions [41]. We validate these results by performing direct numericalsimulations using (6) with outflux of charges through the y boundaries, assuming that theserepresent the charge collectors/electrodes [36]: (cid:0) J uy , J χy , J py , J ny , φ (cid:1) (cid:12)(cid:12) y =0 = (cid:18) , , − D p p ∂φ∂y , , V (cid:19) , (cid:0) J uy , J χy , J py , J ny , φ (cid:1) (cid:12)(cid:12) y = L y = (cid:18) , , , − D n n ∂φ∂y , − V (cid:19) , where V = 0 is a fixed voltage under short circuit conditions, and the fluxes (in their dimen-sionless forms) are J u = D u (1 − u ) ∇ δ E δu , J χ = χ ∇ δ E δχ , J p = D p p ∇ δ E δp , J n = D n n ∇ δ E δn . In the x -direction we employ periodic BC for all fields.At low illumination values, G < G ZZ , we find that the stripes are stable (left inset in Fig. 5).In the region G ZZ < G < G CR , the stripes are unstable only to ZZ, which as can be expecteddevelops in the bulk (middle inset in Fig. 5). The agreement with the linear analysis is excellentand reproduces similar wavenumber k max y = 0 . , as shown by the green curves in the middleinset.In contrast, for G > G CR the primary instability now develops near the y boundaries and is ofa cross-roll mode type (right inset in Fig. 5). Consequently, the charge-flux boundary conditions13reak the degeneracy of the ZZ and the CR modes by enhancing the latter. Nevertheless, theresults are in agreement with the linear analysis (see orange lines near the boundaries) for boththe developed wavenumbers and the onsets (as shown by the dots (symmetric case) and invertedtriangles (asymmetric case)). Consequently, these results indicate that decreasing ξ and thus,pronouncing the mixing energy towards a triple well, shifts the instability onsets to higher G values. The latter in turn, suggests that the OPV will become less susceptible to deformationmodes that enhance morphological degradation, in particular the dangerous CR instability. V. DISCUSSION
Following recent highlights of a three-phase (donor/mixed/acceptor) bulk-heterojunction(BHJ) in organic photovoltaics (OPV) [25, 32, 33, 35], we used a mean-field approach [36] toidentify the role of the intermediate mixed-phase on morphological changes. Under illuminationthe model is driven out of equilibrium so that stripe morphology may arise (Fig. 3). In contrast,under dark conditions the system evolves solely by coarsening [35, 41]. From a mathemati-cal point of view, stripe morphology arises due to a finite wavenumber instability [36] that ispossible only under illumination and whose nature is effected by the order parameter and theexciton/electron/hole fields (see system (6)). We focus on and distinguish between two generictransverse instabilities of donor-acceptor stripes in 2D (distinctly from the formation of stripesby phase separation) with symmetric and asymmetric compositions (as summarized in Fig. 5):the bending (zigzag mode) and the pinching (cross-roll mode). The pinching mode is charac-terized by high mixing energy whereas at low mixing energies bending of the donor/acceptordomains is favored. We emphasize that the time scale separation between the morphological(material) changes and charge dynamics is of several orders of magnitude so that our resultsindicate only the initial trend and not necessarily convergence to a final state, but the furtherevolution, in reality, is extremely slow. Furthermore, the slow time evolution of the materiallowers the sensitivity of the OPV to finite-amplitude perturbations, thus, the effect for exampleof sudden changes in illumination is negligible.Although we limited our analysis to 2D, standard theory shows that the pinching mode mayalso lead to discontinuous and isolated domains in 3D [52–55] and thus, in OPV loss of current tothe electrodes that cause operation failure. This phenomenon resembles the so-called pearling of14ylindrical threads [56–60]. Moreover, according to numerical simulations, relatively large D/Avolumes of BHJ, are more susceptible to transverse instabilities since the intermediate phasedoes not suppress transverse front instabilities that arise due to curvature effects as in bistablesystems [61–64], i.e., in the direction that is parallel to the electrodes. This is consistent withthe diffusion length of about tens of nanometer size of the BHJ [16, 25].Consequently, we showed that the qualitative significance of three-phase BHJ goes beyondinhibition of phase separation [35], as it may have tailoring by demand properties that can becontrolled by the composition of the mixed-phase via donor-acceptor choices: by decreasing themixing energy parameter ξ , the instability onsets are shifted to higher illumination values G .This degree of control is absent or less sequential in two-phase OPV. We believe that our resultsmay assist in the future design of long-lasting OPV, consisting of three-phase BHJ. In a broadercontext, our results should apply to other systems in physicochemical systems that exhibitphase separation [65, 66] and can be driven out of equilibrium, in particular in ion-intercalatedrenewable batteries that depend on reversible phase exchanges in charge/discharge cycling [67–69], such as Li [70–72] and Ni [73–75] based electrodes. ACKNOWLEDGMENTS
The research was done in the framework of the Grand Technion Energy Program (GTEP) andof the BGU Energy Initiative Program, and supported by the Adelis Foundation for renewableenergy research. [1] G. P. Kini, S. J. Jeon, and D. K. Moon, Advanced Materials , 1906175 (2020).[2] H. S. Vogelbaum and G. Sauvé, Synthetic Metals , 107 (2017).[3] A. Bonasera, G. Giuliano, G. Arrabito, and B. Pignataro, Molecules , 2200 (2020).[4] R. Zhou, Z. Jiang, C. Yang, J. Yu, J. Feng, M. A. Adil, D. Deng, W. Zou, J. Zhang, K. Lu, et al. ,Nature Communications , 1 (2019).[5] M. He, F. Qiu, and Z. Lin, Journal of Materials Chemistry , 17039 (2011).[6] J. U. Lee, J. W. Jung, J. W. Jo, and W. H. Jo, Journal of Materials Chemistry , 24265 (2012).[7] N. Grossiord, J. M. Kroon, R. Andriessen, and P. W. M. Blom, Organic Electronics , 432 (2012).
8] B. A. Collins, J. R. Tumbleston, and H. Ade, Journal of Physical Chemistry Letters , 3135 (2011).[9] F. Liu, Y. Gu, J. W. Jung, W. H. Jo, and T. P. Russell, Journal of Polymer Science Part B:Polymer Physics , 1018 (2012).[10] D. R. Kozub, K. Vakhshouri, L. M. Orme, C. Wang, A. Hexemer, and E. D. Gomez, Macro-molecules , 5722 (2011).[11] K. Vakhshouri, D. R. Kozub, C. Wang, A. Salleo, and E. D. Gomez, Physical Review Letters ,026601 (2012).[12] C. J. Schaffer, C. M. Palumbiny, M. A. Niedermeier, C. Jendrzejewski, G. Santoro, S. V. Roth,and P. Müller-Buschbaum, Advanced Materials , 6760 (2013).[13] C. J. Schaffer, C. M. Palumbiny, M. A. Niedermeier, C. Burger, G. Santoro, S. V. Roth, andP. Müller-Buschbaum, Advanced Energy Materials , 1600712 (2016).[14] H. B. Naveed, K. Zhou, and W. Ma, Accounts of Chemical Research , 2904 (2019).[15] W. R. Mateker and M. D. McGehee, Advanced Materials , 1603940 (2017).[16] M. Jørgensen, K. Norrman, S. A. Gevorgyan, T. Tromholt, B. Andreasen, and F. C. Krebs,Advanced Materials , 580 (2012).[17] B. A. Collins, E. Gann, L. Guignard, X. He, C. R. McNeill, and H. Ade, Journal of PhysicalChemistry Letters , 3160 (2010).[18] J. Zhao, A. Swinnen, G. Van Assche, J. Manca, D. Vanderzande, and B. Van Mele, Journal ofPhysical Chemistry B , 1587 (2009).[19] N. D. Treat, M. A. Brady, G. Smith, M. F. Toney, E. J. Kramer, C. J. Hawker, and M. L. Chabinyc,Advanced Energy Materials , 82 (2011).[20] S. Kouijzer, J. J. Michels, M. van den Berg, V. S. Gevaerts, M. Turbiez, M. M. Wienk, and R. A.Janssen, Journal of the American Chemical Society , 12057 (2013).[21] N. D. Treat and M. L. Chabinyc, Annual Review of Physical Chemistry , 59 (2014).[22] I. Cardinaletti, J. Kesters, S. Bertho, B. Conings, F. Piersimoni, J. D’Haen, L. Lutsen, M. Nesladek,B. Van Mele, G. Van Assche, K. Vandewal, A. Salleo, D. Vanderzande, W. Maes, and J. V. Manca,Journal of Photonics for Energy , 040997 (2014).[23] U. Vongsaysy, D. M. Bassani, L. Servant, B. Pavageau, G. Wantz, and H. Aziz, Journal of Photonicsfor Energy , 040998 (2014).[24] K. Zhou, J. Liu, M. Li, X. Yu, R. Xing, and Y. Han, Journal of Physical Chemistry C , 1729 , 4234(2014).[26] O. Reid, J. Malik, G. Latini, S. Dayal, N. Kopidakis, C. Silva, N. Stingelin, and G. Rumbles,Journal of Polymer Science, Part B: Polymer Physics , 27 (2012).[27] J. Razzell-Hollis, W. C. Tsoi, and J.-S. Kim, Journal of Materials Chemistry C , 6235 (2013).[28] J. Bartelt, Z. Beiley, E. Hoke, W. Mateker, J. Douglas, B. Collins, J. Tumbleston, K. Graham,A. Amassian, H. Ade, J. Fréchet, M. Toney, and M. Mcgehee, Advanced Energy Materials , 364(2013).[29] T. Burke and M. McGehee, Advanced Materials , 1923 (2014).[30] P. Müller-Buschbaum, Advanced Materials , 7692 (2014).[31] N. Gasparini, X. Jiao, T. Heumueller, D. Baran, G. J. Matt, S. Fladischer, E. Spiecker, H. Ade,C. J. Brabec, and T. Ameri, Nature Energy , 16118 (2016).[32] K. Zhou, J. Xin, and W. Ma, ACS Energy Letters , 447 (2019).[33] C. Wang, X. Xu, W. Zhang, S. B. Dkhil, X. Meng, X. Liu, O. Margeat, A. Yartsev, W. Ma,J. Ackermann, et al. , Nano Energy , 24 (2017).[34] F. Liu, W. Zhao, J. R. Tumbleston, C. Wang, Y. Gu, D. Wang, A. L. Briseno, H. Ade, and T. P.Russell, Advanced Energy Materials , 1301377 (2014).[35] S. B. Dkhil, M. Pfannmöller, M. I. Saba, M. Gaceur, H. Heidari, C. Videlot-Ackermann, O. Margeat,A. Guerrero, J. Bisquert, G. Garcia-Belmonte, et al. , Advanced Energy Materials , 1601486 (2017).[36] A. Z. Shapira, N. Gavish, and A. Yochelis, EPL (Europhysics Letters) , 38001 (2019).[37] W. Ma, J. R. Tumbleston, M. Wang, E. Gann, F. Huang, and H. Ade, Advanced Energy Materials , 864 (2013).[38] N. Gavish, I. Versano, and A. Yochelis, SIAM Journal on Applied Dynamical Systems , 1946(2017).[39] A. Z. Shapira, H. Uecker, and A. Yochelis, Chaos: An Interdisciplinary Journal of NonlinearScience , 073104 (2020).[40] M. C. Cross and P. C. Hohenberg, Reviews of Modern Physics , 851 (1993).[41] G. Buxton and N. Clarke, Physical Review B , 085207 (2006).[42] H. Uecker and D. Wetzel, SIAM Journal on Applied Dynamical Systems , 94 (2014).
43] T. Dohnal, J. D. Rademacher, H. Uecker, and D. Wetzel, Proceedings of ENOC14 (2014).[44] H. Greenside and W. Coughran Jr, Physical Review A , 398 (1984).[45] H. Greenside and M. Cross, Physical Review A , 2492 (1985).[46] U. Thiele and E. Knobloch, Physics of Fluids , 892 (2003).[47] T. Kolokolnikov, M. J. Ward, and J. Wei, Studies in Applied Mathematics , 35 (2006).[48] T. Kolokolnikov, W. Sun, M. Ward, and J. Wei, SIAM Journal on Applied Dynamical Systems ,313 (2006).[49] J. Burke and E. Knobloch, Chaos , 037102 (2007).[50] J. A. Diez, A. G. González, and L. Kondic, Physics of Fluids , 032104 (2012).[51] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes 3rd edition:The art of scientific computing (Cambridge University Press, 2007).[52] X. Yu and J. T. Liu, Physics of Fluids , 736 (1994).[53] V. V. Kolmychkov, O. S. Mazhorova, Y. P. Popov, P. Bontoux, and M. El Ganaoui, ComptesRendus Mecanique , 739 (2005).[54] E. Fedoseev, V. Kolmychkov, and O. Mazhorova, Progress in Computational Fluid Dynamics, anInternational Journal , 208 (2010).[55] H. Uecker and D. Wetzel, Physica D , 132383 (2020).[56] I. Tsafrir, D. Sagi, T. Arzi, M.-A. Guedeau-Boudeville, V. Frette, D. Kandel, and J. Stavans,Physical Review Letters , 1138 (2001).[57] P. Nelson, T. Powers, and U. Seifert, Physical Review Letters , 3384 (1995).[58] K. P. Sinha, S. Gadkari, and R. M. Thaokar, Soft Matter , 7274 (2013).[59] S. Chaïeb and S. Rica, Physical Review E , 7733 (1998).[60] T. Nguyen, A. Gopal, K. Lee, and T. Witten, Physical Review E , 051930 (2005).[61] R. E. Goldstein, D. J. Muraki, and D. M. Petrich, Physical Review E , 3933 (1996).[62] A. Yochelis, C. Elphick, A. Hagberg, and E. Meron, Physica D , 201 (2004).[63] A. Hagberg, A. Yochelis, H. Yizhaq, C. Elphick, L. Pismen, and E. Meron, Physica D , 186(2006).[64] T. Kolokolnikov and M. Tlidi, Physical Review Letters , 188303 (2007).[65] H. Emmerich, Advances in Physics , 1 (2008).[66] S. DeWitt and K. Thornton, in Computational Materials System Design (Springer, 2018) pp. 67–87.
67] J. L. Kaufman, J. Vinckeviči¯ut˙e, S. Krishna Kolli, J. Gabriel Goiri, and A. Van der Ven, Philo-sophical Transactions of the Royal Society A , 20190020 (2019).[68] A. R. Balakrishna, Y.-M. Chiang, and W. C. Carter, Physical Review Materials , 065404 (2019).[69] A. Van der Ven, Z. Deng, S. Banerjee, and S. P. Ong, Chemical Reviews , 6977 (2020).[70] M. Tang, W. C. Carter, and Y.-M. Chiang, Annual Review of Materials Research , 501 (2010).[71] D. Grazioli, M. Magri, and A. Salvadori, Computational Mechanics , 889 (2016).[72] Y. Zhao, P. Stein, Y. Bai, M. Al-Siraj, Y. Yang, and B.-X. Xu, Journal of Power Sources , 259(2019).[73] G. Briggs and M. Fleischmann, Transactions of the Faraday Society , 2397 (1971).[74] R. Barnard, C. Randell, and F. Tye, Journal of Applied Electrochemistry , 109 (1980).[75] R. Huggins, H. Prinz, M. Wohlfahrt-Mehrens, L. Jörissen, and W. Witschel, Solid State Ionics ,417 (1994).,417 (1994).